【发布时间】:2016-04-05 23:17:22
【问题描述】:
我在 R 中创建了一个新的 Robust HoltWinters 函数(基于 stats::Holt-Winters)方法(根据“使用指数和 Holt-Winters 进行稳健预测 平滑”,作者 Sarah Gelper1,Roland Fried,Christophe Croux。2008 年 9 月 26 日。)为什么?嗯...为什么不呢!但我离题了...
stats::Holt-Winters 方法的核心是一个名为 C_HoltWinters 的 C 代码,我已对其进行了修改以使其更加健壮(见下文)
#include <stdlib.h>
#include <string.h> // memcpy
#include <math.h>
#include <R.h>
#include "ts.h"
void HoltWinters (
double *x, /*as.double(x) */
double *x_adj, /*Adjust time series data, if need be Added*/
int *xl, /*lenx - Length of the current time series*/
double *alpha, /*as.double(max(min(alpha, 1), 0)), */
double *beta, /*as.double(max(min(beta,1), 0)), */
double *gamma, /*as.double(max(min(gamma, 1), 0)), */
double *llamda,/*as.double(max(min(llamda,1),0)), ADDED*/
int *start_time, /*as.integer(start.time), */
int *seasonal, /*as.integer(!+(seasonal == "multiplicative")), */
int *period, /* as.integer(f), */
int *dotrend, /* as.integer(!is.logical(beta) || beta), */
int *doseasonal, /* as.integer(!is.logical(gamma) || gamma), */
double *a, /*l.start - starting values for level*/
double *b, /*b.start - starting values for Trend*/
double *s, /*s.start - starting values for SEasonal*/
double *l, /*t.start - starting values for LLamda ADDED*/
double *k, /* Value for K ADDED*/
double *ck, /*value for ck ADDED*/
/* return values */
double *SSE,
double *level,
double *trend,
double *season
)
{
double res = 0, xhat = 0, stmp = 0, theta = 1, RhoK = 0, phi = 0 ;
int i, i0, s0; /*i is the current t, i0 is the current LESS starting period, and s0 = is the seasonal current LESS Starting period*/
/* copy start values to the beginning of the vectors */
level[0] = *a;
if (*dotrend == 1) trend[0] = *b;
if (*doseasonal == 1) memcpy(season, s, *period * sizeof(double));
for (i = *start_time - 1; i < *xl; i++) {
/* indices for period i */
i0 = i - *start_time + 2;
s0 = i0 + *period - 1;
/* forecast *for* period i */
xhat = level[i0 - 1] + (*dotrend == 1 ? trend[i0 - 1] : 0);
stmp = *doseasonal == 1 ? season[s0 - *period] : (*seasonal != 1);
if (*seasonal == 1)
xhat += stmp;
else
xhat *= stmp;
/* Sum of Squared Errors */
res = x[i] - xhat;
/*adjusting for robustness....Gahds*/
RhoK = (abs(res / theta) <= *k ? *ck * (1 - pow(1 - pow((res / (*k * theta)),2),3)): *ck);
theta = sqrt(*llamda * RhoK * pow(theta,2) + (1 - *llamda) * pow(theta,2));
phi = (abs(res / theta) < *k ? res / theta : ((res / theta) / abs(res / theta) * (*k)));
x_adj[i] = phi * theta + xhat;
res = x_adj[i] - xhat;
*SSE += res * res;
/* estimate of level *in* period i */
if (*seasonal == 1)
level[i0] = *alpha * (x_adj[i] - stmp)
+ (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
else
level[i0] = *alpha * (x_adj[i] / stmp)
+ (1 - *alpha) * (level[i0 - 1] + trend[i0 - 1]);
/* estimate of trend *in* period i */
if (*dotrend == 1)
trend[i0] = *beta * (level[i0] - level[i0 - 1])
+ (1 - *beta) * trend[i0 - 1];
/* estimate of seasonal component *in* period i */
if (*doseasonal == 1) {
if (*seasonal == 1)
season[s0] = *gamma * (x_adj[i] - level[i0])
+ (1 - *gamma) * stmp;
else
season[s0] = *gamma * (x_adj[i] / level[i0])
+ (1 - *gamma) * stmp;
}
}
}
所以我在 windows sigh 中用 R (3.2.2) 编译它:
R CMD SHLIB C_R_HoltWinters.c
gcc -m64 -I"C:/PROGRA~1/R/R-32~1.2/include" -DNDEBUG -I"d:/RCompile/r-compiling/local/local320/include" -O2 -Wall -std=gnu99 -mtune=core2 -c C_R_HoltWinters.c -o C_R_HoltWinters.o
gcc -m64 -shared -s -static-libgcc -o C_R_HoltWinters.dll tmp.def C_R_HoltWinters.o -Ld:/RCompile/r-compiling/local/local320/lib/x64 -Ld:/RCompile/r-compiling/local/local320/lib -LC:/PROGRA~1/R/R-32~1.2/bin/x64 -lR
将其加载到 R 中:
dyn.load('C_R_HoltWinters.dll')
检查是否存在
> getLoadedDLLs()
Filename Dynamic.Lookup
base base FALSE
utils C:/Program Files/RRO/R-3.2.2/library/utils/libs/x64/utils.dll FALSE
methods C:/Program Files/RRO/R-3.2.2/library/methods/libs/x64/methods.dll FALSE
RevoUtilsMath C:/Program Files/RRO/R-3.2.2/library/RevoUtilsMath/libs/x64/RevoUtilsMath.dll TRUE
grDevices C:/Program Files/RRO/R-3.2.2/library/grDevices/libs/x64/grDevices.dll FALSE
graphics C:/Program Files/RRO/R-3.2.2/library/graphics/libs/x64/graphics.dll FALSE
stats C:/Program Files/RRO/R-3.2.2/library/stats/libs/x64/stats.dll FALSE
tools C:/Program Files/RRO/R-3.2.2/library/tools/libs/x64/tools.dll FALSE
internet C:/PROGRA~1/RRO/R-32~1.2/modules/x64/internet.dll TRUE
(embedding) (embedding) FALSE
C_R_HoltWinters C:/scripts/R/C_R_HoltWinters.dll TRUE
啊,是的,就是这样。所以,只是为了大便和咯咯笑,我对其进行了测试:
> is.loaded('C_R_HoltWinters')
[1] FALSE
> is.loaded("C_R_HoltWinters")
[1] FALSE
> is.loaded(C_R_HoltWinters)
Error in is.loaded(C_R_HoltWinters) : object 'C_R_HoltWinters' not found
好吧....它应该在那里,但它不是。也许它知道我不知道的东西,所以我尝试运行它:
> .C("C_R_HoltWinters", blahblahblah)
Error in .C("C_R_HoltWinters") :
C symbol name "C_R_HoltWinters" not in load table
> .Call("C_R_HoltWinters", blahblahblah)
Error in .Call("C_R_HoltWinters") :
C symbol name "C_R_HoltWinters" not in load table
但是当我加载一个名为 foo 的不同 c 代码并运行它时,它运行良好。
为什么 R 不能引用 C_R_HoltWinters.dll?如果我把它放在一个包裹里,这也会坏吗?
谢谢
【问题讨论】:
-
没什么,执行得很好
-
这真的很奇怪。你能在
getLoadedDLLs()的输出中找到这个 HoltWinters 吗?另一方面,如果您正在考虑在包中使用此代码,我强烈建议查看 Rcpp,因为它更容易导航外部调用在 R 中的混乱。 -
嗯。我在这里编译了你的 C 代码,删除了
#include "ts.h"行。它成功dyn.loads,出现在getLoadedDLLs(),但是当我检查它是否is.loaded时,它返回FALSE。奇怪。 -
PS,要运行它,您应该使用函数的名称,而不是文件。所以应该是
.C("HoltWinters", bla bla bla),这显然是在为我运行。 -
是的,我可以在 getLoadedDlls() 中找到 C_R_HoltWinters,其值为 Dynamic.Lookup = True
标签: c r statistics