在编程中，有两样可以减少代码重复的工具，一是循环，一是函数。本文用到的循环主要是for循环和while循环。python中可以用def定义函数，sas中可以用宏写函数，R中可以用function定义函数。
任务：生成2*2的随机数矩阵，进行矩阵运算，模拟100次，输出模拟试验的计数值与结果，同时保证每次的计算结果中没有无效值。
思路：
NO.1 将矩阵运算封装到一个函数内，即定义一个矩阵运算的函数，这样输入原始值后，就可以输出一次模拟值；
NO.2 建一个空的数据框，设置种子，将函数放入循环中，在循环中调用函数，循环多少次就会输出多少次模拟值，并且将循环结果写到数据框中，就可以得到一个模拟100次矩阵计算的数据集。
函数基础知识：
1、R语言笔记7：functions——编写函数所需的基础知识 - 百度文库 (baidu.com)
2、2 IT - UNIT-1 Start Learning R.pdf (gvpcew.ac.in)
3、140.776 Statistical Computing课程： 1）StatisticsFunc2.pdf (jhsph.edu)；2）StatFunc.pdf (jhsph.edu)
程序与分析
#生成2*2的随机数矩阵：Q<-matrix(sample(100:200,4,replace=TRUE),ncol=2)
#函数计算中用到的条件
IOratio<-c(0.7,0.75)
l<-c(0.01,0.02)
k<-c(0.4,0.5)
L<-c(0.2,0.1)
K<-c(0.2,0.1)
I<-diag(2)
#定义矩阵计算的函数，此处要注意的函数的初始条件，可以是一个数，可以是数组，可以是字符串，也可以是矩阵，以下代表的就是给定初始矩阵矩阵M，计算由其生成的相关值。
#函数里的print和return有什么区别呢？假如没有print也没有return,则输出函数中的最后一个结果；print代表在输出台中打印部分结果，不影响运行函数的后续部分。return表示返回结果，即可以控制输出函数的哪一部分。
fun<- function(M) {
MI <-colSums(M[,])
MI <-data.matrix(MI)
X<-MI/IOratio
V<-X-MI
ln_V<-log(V)
L_V<-L*X/V
K_V<-K*X/V
l_s<-l*X
k_s<-k*X
ln_l_s<-log(l_s)
ln_k_s<-log(k_s)
V_T<-colSums(V)
W_V<-t(V/V_T)
TFP_s_V<-ln_V-L_V*ln_l_s-K_V*ln_k_s
APG_V<-W_V %*% TFP_s_V
id<-i
print(id)
test<- cbind.data.frame(id, APG_V, TFP_s_V[1,1], TFP_s_V[2,1])
return(test)
}
#建立一个空的数据框，存储循环得到的结果。以后可以讲解利用列表存储更复杂的循环结果。
hh<-data.frame()
#利用for循环实现进行100次有效的随机试验。for(i in 1:100)的解读：进行100次试验；在最后写成的数据框里，每一行就是一次试验，因为我们在之前的函数中，引入了id，所以输出的数据集中，id记录了是第几次试验。
#利用while 判断每次实验的结果符不符合条件。while(sum(is.na(test1))>0)的解读：只要sum(is.na(test1))>0，即test1里有一个或多个缺失值，就会继续while循环。如果符合条件，则不储存结果，继续运行函数，直到不符合条件，才认为是一次有效的实验，并存储实验结果。注：sum(is.na(test1))>0是判断一行里有没有缺失值，is.na(XX$XXX)=TRUE是判断某行的某单个变量是不是缺失值。
#由于在循环之前建造了一个空的数据集，所以我们每进行一次试验，就会给这个数据集填上新的一行。每多进行一次试验，就会在之前的结果后面多记录一次新实验的结果。
#set.seed():一个特定的种子可以产生一个特定的伪随机序列。引入种子可以让模拟能可重复出现，因为很多时候我们需要取随机数，但这段代码再跑一次的时候，结果就不一样了，如果需要重复出现同样的模拟结果的话，就可以用set.seed()。
for(i in 1:100)
{set.seed(i)
test1=as.data.frame(fun(matrix(sample(100:200,4,replace=TRUE),ncol=2)))
while(sum(is.na(test1))>0){
test1=as.data.frame(fun(matrix(sample(100:200,4,replace=TRUE),ncol=2)))
}
hh<- rbind(hh,test1)
}
网友评论