发帖
雷达卡
第三节 R的最适化运算
R语言内建了非线性最适化的运算功能，这在模型参数的校准上非常有用。R提供了nlm()与optim()这两个函数，让我们进行未受限的最适化。使用者可以选择使用BFGS、共轭梯度、退火模拟等方法，进行计算。
一、简单最适化案例
我们先介绍nlm()的使用，nlm()需要传入一个目标函数，起始向量，以及目标函数中的参数值。令目标函数如下，
其中，
为搜寻的变量值，A、B为参数。我们希望在给定A、B参数值下，找到使目标函数值最小化的变量值。
......................................................................(4.3.2)
> g<- function(x, A, B) {
+ out <- sin(x[1]) - sin(x[2]-A) + x[3]^2+ B
+ return(out) }
> results <- nlm(g, c(1,2,3), A=4,B=2)
> results
$minimum
[1] 6.497025e-13
$estimate
[1] -1.570797e+00 -7.123895e-01-4.990548e-07
$gradient
[1] 5.300922e-08 -1.465494e-08 1.776357e-09
$code
[1] 1
$iterations
[1] 12
程序行表4.3.1
上面程序先建立目标函数g，呼叫nlm()时，需传入目标函数g，猜测的初始值向量c(1,2, 3)，以及给定的A、B参数值。极小值为6.497025*10-13，估计的最适解为(-1.570797,-0.7123895, -4.990548*10-7)。
二、受限最适化与惩罚函数法
对于受限函数的最适化，一个简单的做法就是在未受限的目标函数，加上惩罚函数。所谓的惩罚函数就是针对超出限制范围的变量，给予目标函数偏离极大的数值。举例而言，如果上例中我们希望x2>=0，我们可以设定目标函数为，
如果x2<0，则惩罚函数将产生一个很大的正数，此很大的正数将不利于目标函数的极小化。
> h <- function(x, A, B) {
+ if(x[2]>=0) res <- g(x, A, B) else res <- g(x, A, B)-1000000*x[2]
+ return(res) }
> results <- nlm(h, c(1,2,3), A=4,B=2)
> results
$minimum
[1] 0.3501565
$estimate
[1] -1.106889e+00 7.153626e-06 3.557386e-02
$gradient
[1] 0.44744613 0.65364941 0.07114872
$iterations
[1] 26
> g(c(1, 1, 1), A=4, B=2)
[1] 3.982591
程序行表4.3.2
上面程序中，我们要求第二个变量的解要大于0，的最适估计值为0.3501565。估计的最适解为(-1.106889,7.153626*10-6, 3.557386*10-2)。
第四节 R的随机数与统计分配函数
R已内建多个分配的随机数生成器，这对我们进行模拟非常便利。下面概述这些随机数生成器使用的方式。
一、设定随机数种子
在使用随机数生成器时，首先要设定随机数种子，我们可以使用下面指令，
set.seed(seed, kind = NULL, normal.kind =NULL)
第一个参数seed为一整数，是我们要设定的随机数种子。
第二个参数kind为一字符串或NULL，是我们要设定的随机数生成器种类。可以选择的种类字符串有，"Wichmann-Hill"、"Marsaglia-Multicarry"、"Super-Duper"、"Mersenne-Twister"、"Knuth-TAOCP-2002"、"Knuth-TAOCP"、"L'Ecuyer-CMRG"、"user-supplied"，若没有设定，预设为"Mersenne-Twister"。
第三个参数normal.kind为一字符串或NULL，是我们要设定的产生常态分配随机数的方法。可以选择的种类字符串有，"Kinderman-Ramage"、"Buggy Kinderman-Ramage"、"Ahrens-Dieter"、"Box-Muller"、"Inversion"、"user-supplied"，若没有设定，预设为"Inversion"。
使用范例
set.seed(1234);
二、均等分配随机数
Usage
runif(n,min = 0, max = 1)
Arguments
| x, q | vector of quantiles. |
| p | vector of probabilities. |
| n | number of observations. If length(n) > 1, the length is taken to be the number required. |
| min, max | lower and upper limits of the distribution. Must be finite. |
| log, log.p | logical; if TRUE, probabilities p are given as log(p). |
| lower.tail | logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x]. |
对应的统计分配函数
dunif(x,min = 0, max = 1, log = FALSE)：机率密度函数(pdf)
punif(q,min = 0, max = 1, lower.tail = TRUE, log.p = FALSE)：累积机率密度函数(cdf)
使用范例
set.seed(1234);
OBS = 1024;
UnifRV = numeric(OBS);
UnifRV = runif(OBS);
print(UnifRV);
三、常态分配随机数
Usage
rnorm(n,mean = 0, sd = 1)
dnorm(x,mean = 0, sd = 1, log = FALSE)
pnorm(q,mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
Arguments
| x, q | vector of quantiles. |
| P | vector of probabilities. |
| n | number of observations. If length(n) > 1, the length is taken to be the number required. |
| mean | vector of means. |
| sd | vector of standard deviations. |
| Log, log.p | logical; if TRUE, probabilities p are given as log(p). |
| lower.tail | logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x]. |
使用范例
set.seed(1234);
OBS = 1024;
NormRV = numeric(OBS);
NormRV = rnorm(OBS);
print(NormRV);
四、Poisson分配随机数
Usage
rpois(n,lambda)
dpois(x,lambda, log = FALSE)
ppois(q,lambda, lower.tail = TRUE, log.p = FALSE)
Arguments
| x | vector of (non-negative integer) quantiles. |
| q | vector of quantiles. |
| p | vector of probabilities. |
| n | number of random values to return. |
| lambda | vector of (non-negative) means. |
| log, log.p | logical; if TRUE, probabilities p are given as log(p). |
| lower.tail | logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x]. |
使用范例
set.seed(1234);
OBS = 1024;
PoisRV = numeric(OBS);
PoisRV = rpois(OBS);
print(PoisRV);
京公网安备 11010802022788号
