重复运行一个r程序多次

jiachunyang1988 发表于 2014-12-2 20:33
我晓得你的意思了，可是我的程序是分为两个部分，第一部分估计得到一个参数的值，在得到第一个估计值的基 ...

不好意思啊，麻烦你帮我看看
###生成数据

error=c()
sigma=1
error=rnorm(20, mean=0, sd=sigma)
error

beta= c()
r0=1
r1=0.4
alpha=1
T=20

beta[1]=error[1]
for (t in 2:length(error)) {
beta[t]=r0+r1*beta[t-1]+error[t]
}
beta=ts(beta)
beta  


X=matrix(runif(600,0,1),nrow=20)
X

Z=matrix( )
Z=diag(beta)%*%X
Z  

A=matrix(c(alpha),nrow=20,ncol=30)+Z
A

P=exp(A)/(1+exp(A))
P

W=matrix(runif(600,0,1),20,30)  
Y=1*(W<P)
Y   
#####################################################################

####给alpha和beta一个初值

a=matrix(0,20,2)
   for(t in 1:T){
   b=glm(Y[t,]~X[t,],family=binomial)
   a[t,]=b$coef
}
a
alpha1=mean(a[,1])
alpha1
beta1=a[,2]
beta1   
###############################################################################

###牛顿法得到非线性方程组的解 r0,r1,sigma

Newtons=function(fun,phi, ep=1e-5, it_max=100){
      
   

   index=0; k=1

   while (k<=it_max){
     
      phi1=phi; obj=fun(phi);

      phi=phi-solve(obj$J, obj$f);

      norm=sqrt((phi-phi1)%*%(phi-phi1))

      if(norm<ep){

         index=1; break

      }
      k<-k+1

   }

   obj=fun(phi);

   list(root=phi, it=k, index=index, FunVal=obj$f)

}

t=2:T

funs=function(phi){
   
   f=c((1+phi[3])*(beta1[1]-phi[1]/(1-phi[3]))/phi[2]^2+(sum(beta1[t]-phi[1]-phi[3]*beta1[t-1]))/phi[2]^2,
   
      -T/phi[2]+(sum((beta1[t]-phi[1]-phi[3]*beta1[t-1])^2))/phi[2]^3+((1-phi[3]^2)*(beta1[1]-phi[1]/(1-phi[3]))^2)/phi[2]^3,
     
      - phi[3]/(1-phi[3]^2)+phi[3]*((beta1[1]-phi[1]/(1-phi[3]))^2)/phi[2]^2+
      
       (sum((beta1[t]-phi[1]-phi[3]*beta1[t-1])*beta1[t-1]))/phi[2]^2+phi[1]*(1+phi[3])*(beta1[1]-phi[1]/(1-phi[3]))/(phi[2]^2*(1-phi[3])))


     J=matrix(c(-(1+phi[3])/(phi[2]^2*(1-phi[3]))-(T-1)/phi[2]^2,
           
            -2*(1+phi[3])*(beta1[1]-phi[1]/(1-phi[3]))/phi[2]^3-2*(sum(beta1[t]-phi[1]-phi[3]*beta1[t-1]))/phi[2]^3,
           
            (beta1[1]*(1-phi[3])^2-2*phi[1])/(phi[2]^2*(1-phi[3])^2)-(sum(beta1[t-1]))/phi[2]^2,
           
            (-2/phi[2]^3)*((sum(beta1[t]-phi[1]-phi[3]*beta1[t-1]))+(1+phi[3])*(beta1[1]-phi[1]/(1-phi[3]))),
            
            T/phi[2]^2-(3/phi[2]^4)*((sum((beta1[t]-phi[1]-phi[3]*beta1[t-1])^2))+(1-phi[3]^2)*(beta1[1]-phi[1]/(1-phi[3]))^2),           
            
            -2*(sum((beta1[t]-phi[1]-phi[3]*beta1[t-1])*beta1[t-1]))/phi[2]^3-2*phi[3]*(beta1[1]-(phi[1]/(1-phi[3]))^2)/phi[2]^3-
            
            2*phi[1]*(1+phi[3])*(beta1[1]-(phi[1]/(1-phi[3])))/(phi[2]^3*(1-phi[3])),
            
            (-2*phi[3]*(beta1[1]-phi[1]/(1-phi[3])))/(phi[2]^2*(1-phi[3]))-(sum(beta1[t-1]))/phi[2]^2+
               
                 (1+phi[3])*(beta1[1]-2*phi[1]/(1-phi[3]))/(phi[2]^2*(1-phi[3])),           
         
             (-2*phi[3]/phi[2]^3)*(beta1[1]-(phi[1]/(1-phi[3]))^2)-2*(sum((beta1[t]-phi[1]-phi[3]*beta1[t-1])*beta1[t-1]))/phi[2]^3-

             (2*phi[1]*(1+phi[3])*(beta1[1]-(phi[1]/(1-phi[3]))))/(phi[2]^3*(1-phi[3])),
               
            -(1+phi[3]^2)/(1-phi[3]^2)^2+(beta1[1]-phi[1]/(1-phi[3]))^2/phi[2]^2-2*phi[1]*phi[3]*(beta1[1]-phi[1]/(1-phi[3]))/(phi[2]^2*(1-phi[3])^2)-
                 
                 (sum(beta1[t-1]^2))/phi[2]^2+(2*phi[1]*beta1[1]-phi[1]^2-((2*phi[1]^2*(1+phi[3]))/(1-phi[3])))/(phi[2]^2*(1-phi[3])^2)),3,3,byrow=T)
   

      list(f=f,J=J)
   }

Newtons(funs,c(1,1,0.4))

能把你的程序贴好点吗。。。看得实在是头痛。。。

roota<-matrix(30,30,3)
inta<-matrix(30,30,1)
indexa<-matrix(30,30,1)
FunVala<-matrix(30,30,3)
for(fi in 1:30)
{
source("test.r")$value->result
result$root->roota[fi,]
result$it->inta[fi,]
result$index->indexa[fi,]
result$FunVal->FunVala[fi,]
}
但是会有警告。。。你这不只四个值吧。。。是四个参数，每个参数的值不一样。。。index好像都是1。。。
resultreal<-cbind(roota,inta,indexa,FunVala)
colnames(resultreal)<-c(paste("root",1:3,sep="")"it","index",paste("FunVal",1:3,sep=""))
resultreal
              root1        root2        root3 it index       FunVal1
[1,]  2.932705e+00     4.443398 -0.044705572 11     1  1.942890e-16
[2,]  1.763472e+00     2.566635  0.131168839  9     1  2.220446e-16
[3,]  1.446558e+00     2.795726  0.072901617 11     1  7.993606e-14
[4,]  3.105975e+00     4.483115  0.111386883 11     1  6.938894e-16
[5,]  2.797802e+00     4.729219  0.161051839 11     1  3.042011e-14
[6,]  2.474712e+00     1.854530 -0.068048059  8     1 -6.800116e-16
[7,]  5.876404e+00     8.698359 -0.218396051 14     1  1.526557e-16
[8,]  3.394085e+00     5.108472 -0.170423350 11     1  6.861178e-14
[9,]  1.570930e+00     2.308090  0.221471928  9     1  6.983303e-14
[10,]  3.167952e+00     7.189025  0.010838923 12     1  2.484124e-15
[11,]  2.930424e+00     3.320860  0.115998544 10     1  7.771561e-16
[12,]  1.380893e+00     2.895080  0.335350108  9     1  3.591571e-14
[13,]  1.798078e+00     1.742767 -0.252371537  9     1  1.809664e-14
[14,]  1.393069e+00     4.084480  0.228269922 11     1  1.360023e-14
[15,]  1.613741e+00     3.408410  0.282311045 10     1 -1.804112e-16
[16,]  4.955300e+01   195.527915 -0.050350331 24     1 -8.673617e-19
[17,]  2.536230e+00     5.259663  0.111575643 11     1  3.143319e-15
[18,]  1.653942e+02   675.067117 -0.050933049 29     1  1.084202e-18
[19,]  3.390971e+00     5.246574  0.226605727 12     1 -1.249001e-16
[20,]  3.618624e+00     3.613593 -0.257438210 11     1 -8.326673e-17
[21,] -5.753179e-03     2.416770  0.362227406  9     1 -4.996004e-16
[22,]  2.961887e+00     2.267074 -0.281913297 10     1  1.465494e-14
[23,]  2.782912e+00     2.137948 -0.177550800  9     1  1.029177e-13
[24,]  1.994353e+00     3.286657  0.002816165  9     1  5.447864e-13
[25,]  1.812346e+00     2.191288  0.008411907  8     1  2.220446e-16
[26,]  1.898667e+00     3.705085 -0.083605849 10     1  2.386980e-14
[27,]  4.262193e-01     1.901698  0.438736775  8     1  5.329071e-15
[28,]  5.100577e+03 20973.343390 -0.052902398 41     1 -2.202286e-20
[29,]  3.388581e+01   127.349666 -0.039546251 23     1 -2.168404e-19
[30,]  2.954689e+00     5.778763 -0.062099543 12     1 -1.387779e-17
           FunVal2       FunVal3
[1,] 1.527827e-10  2.085973e-10
[2,] 8.094432e-09 -2.761690e-08
[3,] 4.388331e-07 -6.477192e-07
[4,] 1.003654e-08 -5.862582e-08
[5,] 1.103887e-07 -6.031496e-07
[6,] 1.185892e-09  2.608705e-09
[7,] 7.098080e-10  8.094551e-09
[8,] 8.393417e-08  1.147064e-07
[9,] 4.203914e-07 -5.463229e-07
[10,] 2.480749e-10 -2.823004e-09
[11,] 5.461260e-08 -2.170808e-07
[12,] 1.697435e-07 -4.729970e-07
[13,] 5.115933e-07  7.090607e-07
[14,] 7.828493e-08 -2.025988e-07
[15,] 5.542269e-10 -2.131502e-09
[16,] 6.956582e-13  2.140046e-11
[17,] 5.927087e-09 -7.292370e-09
[18,] 5.340644e-16  2.551084e-13
[19,] 2.963673e-09 -1.519648e-08
[20,] 4.210173e-09  1.527459e-08
[21,] 9.903663e-09 -3.222366e-08
[22,] 5.359805e-07  1.151900e-06
[23,] 6.872364e-07  5.405425e-07
[24,] 4.830573e-09 -3.470433e-08
[25,] 1.272426e-09 -8.095314e-10
[26,] 7.415512e-08  1.332113e-07
[27,] 2.646835e-08 -6.673648e-08
[28,] 9.949249e-18  2.004785e-13
[29,] 1.783610e-14  5.994094e-13
[30,] 6.892791e-11  1.572562e-10
这样解决你的问题了？

大致就是这样。。。不过我没看懂你的程序是想干嘛。。。

曲散人终 发表于 2014-12-2 22:56
rootaindexa[fi,]
result$FunVal->FunVala[fi,]
}

谢谢，能加你扣扣吗？

jiachunyang1988 发表于 2014-12-3 10:47
谢谢，能加你扣扣吗？

我明天考试了，840212114，没什么大问题的话就明天考完再说吧。。。你都能写出这么繁琐的程序，怎么会解决不了这么一个小问题。。。

那就将replicate中的表达式换成你自己编的函数，在函数中设好需要的参数。
你难道没用过lapply吗，跟它一样的道理。
还需要我讲得更清楚点吗？

曲散人终 发表于 2014-12-1 22:05
setwd()设置工作目录，你可以用getwd()来查看当前工作目录，然后setwd()到你r程序所在目录，接着用这个循环 ...

大神 能否解答下这个source("youfile.r")$value中 value 指的是代码最后输出的值，还是其他值呢？

zhangbo19h 发表于 2015-7-24 10:39
大神 能否解答下这个source("youfile.r")$value中 value 指的是代码最后输出的值，还是其他值呢？

你可以用source（“yourfile.r")->a来看看具体值，不同程序不一样。。。

曲散人终 发表于 2015-7-24 11:18
你可以用source（“yourfile.r")->a来看看具体值，不同程序不一样。。。

我知道不同的程序不一样。source（）是读取R代码的，你在给出的解答程序里，是用$value来引用了一个值，我只是想问一下，你这个value具体指的是哪个值？是程序里的最终结果，还是什么。谢谢你。