2011-08-29 84 views
2

我有以下數據集(CEU):R,plyr,具有複雜的功能

group x  y 
1  -23  100 
1  -0.90 69.62 
1  -0.90 72.03 
2  -23  100 
2  0.69 48.01 
2  0.69 45.63 

對於組中的每個值,我想申請下面指出的x和y值的每個子集的功能。然後,我想將所有結果合併,並將它們寫入一個表中以導出。

我不確定如何應用plyr函數來做到這一點...如果這確實是正確的行爲。

x<-c(-23.0000,-0.9031,-0.9031) 
y<-c(100,85.72,86.65) 

par<-c(16.88,100.28,-.75,4.129) 

dcrit<-function(d) { 
    sumsq<-0 
    for (i in 1:length(x)){ 
     sumsq<-sumsq+ (y[i]-(par[1]+(par[2]-par[1])/(1+10^((x[i]-par[3])*d))))^2  
    } 
    sumsq 
} 

S<-function(par) { 
    a<-par[1] 
    b<-par[2] 
    c<-par[3] 
    d<-par[4] 
    sumsq<-0 
    for (i in 1:length(x)){ 
     sumsq<-sumsq+ (y[i]-(a+(b-a)/(1+10^((x[i]-c)*d))))^2  
    } 
    sumsq 
} 
optim(par,S) 

CEU <- read.csv(file="C:/files/CEU.csv",head=TRUE,sep=",") 
CEU 

data <- ddply(CEU,.(group),function(xy) 
{ 
par[1]<-min(y) 
par[2]<-100 
par[3]<-x[[which.min(abs(y-50))]] 
par[4]<-optimize(dcrit,interval=c(-100,100))$minimum 

o<-optim(par,S) 
par<-o$par 

a<-par[1]; 
b<-par[2]; 
c<-par[3]; 
d<-par[4]; 

k<-(b-a)/(20-a)-1 
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA 
ec20 

z<-(b-a)/(50-a)-1 
if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA 
ec50 

j<-(b-a)/(80-a)-1 
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA 
ec80 

data.frame(ec20, ec50, ec80) 

}) 

data 

的代碼運行沒有錯誤,但僅允許在原始x和y值被設置:

x<-c(-23.0000,-0.9031,-0.9031) 
y<-c(100,85.72,86.65) 

在數據集中的CEU x和y值不使用ddply。它們不會以迭代方式替換原始x和y,因爲它們與組值相關。數據具有適當的組數,ec20/ec50/ec80值準確,但僅適用於原始x和y。

> data 
    group  ec20  ec50  ec80 
1  1 -0.3652977 -0.6843279 -0.8530892 
2  2 -0.3652977 -0.6843279 -0.8530892 
3  3 -0.3652977 -0.6843279 -0.8530892 
4  4 -0.3652977 -0.6843279 -0.8530892 
5  5 -0.3652977 -0.6843279 -0.8530892 
+0

Optimize獲取「f」(函數)的第一個參數和「interval」(範圍)的第二個參數。但是你似乎正在發送一個未定義的函數,'dcrit',然後對結果做些什麼,但是,因爲「S」只出現在你的代碼中。 –

+0

使用完整的代碼編輯原始帖子。謝謝! – Sash

回答

3

它看起來像你有權利,你只需要產生輸出。

我猜這是你的輸出?

k<-(b-a)/(20-a)-1 
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA 
ec20 

z<-(b-a)/(50-a)-1 
if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA 
ec50 

j<-(b-a)/(80-a)-1 
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA 
ec80 

把它們放入一個data.frame在函數的末尾:

... 
    data.frame(ec20, ec50, ec80) 
} 

現在你會得到與他們的data.frame,有三列ec20ec50ec80


你的問題與優化:我認爲問題在於

R中
par[3]<-x[which.min(abs(y-50))] 

[不規整標 - 它得到一個切片 - 在這種情況下data.frame列。該行將par從數字向量變爲list。添加更多括號:

par[3]<-x[[which.min(abs(y-50))]] 
+1

如果這是正確的,你確實比我更好的讀者! ;) – joran

+2

他很勇敢。忽略諸如定義函數和測試等無關緊要的細節。切入追逐。讓我們看看......如果我給他的回答加上一個加號,並且你的評論加一個-1,(不完整)的提問者......我要出多少? –

+0

感謝您的評論,我意識到目前還不清楚我最初想用什麼樣的優化功能。我關心的不是我是否正確使用了優化/直流功能,我知道這些工作。我只是想用ddply做一個迭代的方式。唉,我得到以下錯誤在優化(參數,S)錯誤:(列表)對象不能被強制輸入'雙' – Sash