[專題大討論]R or Python, Who Is Better for Statistical Analysis?

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

R or Python, Who Is Better for Statistical Analysis?

1. class NLS:
   
2.     ''' This provides a wrapper for scipy.optimize.leastsq to get the relevant output for nonlinear least squares. Although scipy provides curve_fit for that reason, curve_fit only returns parameter estimates and covariances. This wrapper returns numerous statistics and diagnostics'''
   
3. 
   
4.     import numpy as np
   
5.     import scipy.stats as spst
   
6.     from scipy.optimize import leastsq
   
7. 
   
8.     def __init__(self, func, p0, xdata, ydata):
   
9.         # Check the data
   
10.         if len(xdata) != len(ydata):
    
11.             msg = 'The number of observations does not match the number of rows for the predictors'
    
12.             raise ValueError(msg)
    
13. 
    
14.         # Check parameter estimates
    
15.         if type(p0) != dict:
    
16.             msg = "Initial parameter estimates (p0) must be a dictionry of form p0={'a':1, 'b':2, etc}"
    
17.             raise ValueError(msg)
    
18. 
    
19.         self.func = func
    
20.         self.inits = p0.values()
    
21.         self.xdata = xdata
    
22.         self.ydata = ydata
    
23.         self.nobs = len( ydata )
    
24.         self.nparm= len( self.inits )
    
25. 
    
26.         self.parmNames = p0.keys()
    
27. 
    
28.         for i in range( len(self.parmNames) ):
    
29.             if len(self.parmNames[i]) > 5:
    
30.                 self.parmNames[i] = self.parmNames[i][0:4]
    
31. 
    
32.         # Run the model
    
33.         self.mod1 = leastsq(self.func, self.inits, args = (self.xdata, self.ydata), full_output=1)
    
34. 
    
35.         # Get the parameters
    
36.         self.parmEsts = np.round( self.mod1[0], 4 )
    
37. 
    
38.         # Get the Error variance and standard deviation
    
39.         self.RSS = np.sum( self.mod1[2]['fvec']**2 )
    
40.         self.df = self.nobs - self.nparm
    
41.         self.MSE = self.RSS / self.df
    
42.         self.RMSE = np.sqrt( self.MSE )
    
43. 
    
44.         # Get the covariance matrix
    
45.         self.cov = self.MSE * self.mod1[1]
    
46. 
    
47.         # Get parameter standard errors
    
48.         self.parmSE = np.diag( np.sqrt( self.cov ) )
    
49. 
    
50.         # Calculate the t-values
    
51.         self.tvals = self.parmEsts/self.parmSE
    
52. 
    
53.         # Get p-values
    
54.         self.pvals = (1 - spst.t.cdf( np.abs(self.tvals), self.df))*2
    
55. 
    
56.         # Get biased variance (MLE) and calculate log-likehood
    
57.         self.s2b = self.RSS / self.nobs
    
58.         self.logLik = -self.nobs/2 * np.log(2*np.pi) - self.nobs/2 * np.log(self.s2b) - 1/(2*self.s2b) * self.RSS
    
59. 
    
60.         del(self.mod1)
    
61.         del(self.s2b)
    
62.         del(self.inits)
    
63. 
    
64.     # Get AIC. Add 1 to the df to account for estimation of standard error
    
65.     def AIC(self, k=2):
    
66.     return -2*self.logLik + k*(self.nparm + 1)
    
67. 
    
68.     del(np)
    
69.     del(leastsq)
    
70. 
    
71.     # Print the summary
    
72.     def summary(self):
    
73.         print
    
74.         print 'Non-linear least squares'
    
75.         print 'Model: ' + self.func.func_name
    
76.         print 'Parameters:'
    
77.         print " Estimate Std. Error t-value P(>|t|)"
    
78.         for i in range( len(self.parmNames) ):
    
79.                 print "% -5s % 5.4f % 5.4f % 5.4f % 5.4f" % tuple( [self.parmNames[i], self.parmEsts[i], self.parmSE[i], self.tvals[i], self.pvals[i]] )
    
80.         print
    
81.         print 'Residual Standard Error: % 5.4f' % self.RMSE
    
82.         print 'Df: %i' % self.df

复制代码

1. import pandas as pd
   
2. 
   
3. respData = pd.read_csv(my datafile here)
   
4. respData['respDaily'] = respData['C_Resp_Mass'] * 24
   
5. 
   
6. # Create the Arrhenius temperature
   
7. 
   
8. respData['Ar'] = -1 / (8.617 * 10**-5 * (respData['Temp']+273))
   
9. 
   
10. # First, define the likelihood null model
    
11. def nullMod(params, mass, yObs):
    
12.     a = params[0]
    
13.     c = params[1]
    
14. 
    
15.     yHat = a*mass**c
    
16.     err = yObs - yHat
    
17.     return(err)
    
18. 
    
19. p0 = {'a':1, 'b':1}
    
20. 
    
21. tMod = NLS(nullMod, p0, respData['UrchinMass'], respData['respDaily'] )
    
22. 
    
23. tMod.summary()
    
24. 
    
25. tMod.AIC()
    
26. 
    
27. tMod.logLik

复制代码

1. The above NLS code is the first that I know if in Python to do statistical NLS (not just curve fitting) and get the output I needed. But, it wasn’t easy, it took me about a week of my off (and on) hours. The following R code to do the same takes me maybe 15 minutes:
   
2. 
   
3. # Import the data
   
4. respiration <- read.csv( my file)
   
5. # Standardize to 24 h
   
6. respiration$C_Resp_Day <- respiration$C_Resp_Mass*24
   
7. # Create the Arrhenius temperature
   
8. respiration$Ar<--1/(8.617*10^-5*(respiration$Temp+273))
   
9. 
   
10. Null.mod <- nls(C_Resp_Day~a*UrchinMass^c, data=respiration, start=list(a=exp(6),c=0.25))
    
11. 
    
12. summary(Null.mod)
    
13. AIC(Null.mod)
    
14. logLik(Null.mod)

复制代码

R is just so much easier for real data analysis and has more functions. Python can do these things, but the modules are scattered (there’s at least three separate modules to fit curves that people have written to do different things) and don’t always give the needed output. The NLS example using Python above is, more or less, a replica of R’s NLS output and the missing pieces can be easily added.

R is currently head-and-shoulders above Python for data analysis, but we remain convinced that Python CAN catch up, easily and quickly. It is entirely possible to do statistical analysis in Python if you want to spend the time coding the analyses yourself. By and large, R is the way to go (for now). We would not be surprised that, if 10 years or so, Python’s data analyses libraries coalesce and evolve into something that can rival R.

For those of you interested in trying this out, the data can be downloaded here.

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Statistical statistica statistic Analysis Analysi provides relevant returns import reason

Lisrelchen 发表于 2014-7-1 02:22
Nonlinear Regression using Python

R is just so much easier for real data analysis and has more fu ...

齐头并进

研究研究

如果只是STAT的话那当然是R了，我不知道如果算像logistic regression这样的模型Python有什么有效地办法。如果是Data Mining，两个差不多吧。但是如果是Numerical Analysis 我觉得还是Python方便一些

charlieq 发表于 2014-7-1 04:02
如果只是STAT的话那当然是R了，我不知道如果算像logistic regression这样的模型Python有什么有效地办法。如 ...

那个，logistic regression不是data mining的分类算法之一么。。。。一直搞不懂

不太了解Python，但是觉得R软件不错。貌似也是现在用户增长量增加最快的。R功能较齐全，还能根据自己研究添加程序包，网上也能下载资源，感觉比较强大吧。

不太了解Python，但是觉得R软件不错。貌似也是现在用户增长量增加最快的。R功能较齐全，还能根据自己研究添加程序包，网上也能下载资源，感觉比较强大吧。

complicated 发表于 2014-7-1 09:32
那个，logistic regression不是data mining的分类算法之一么。。。。一直搞不懂

logit regression  是regression的，属于统计方法，Data Mining又是另一回事，两者有思想交叉的地方，但是严格来讲是两个领域。Data Mining主要是大数据处理

charlieq 发表于 2014-7-1 22:21
logit regression  是regression的，属于统计方法，Data Mining又是另一回事，两者有思想交叉的地方，但是 ...

OK~

数据一多，总感觉python好用些。