/* Authors Martin Schlather, schlather@math.uni-mannheim.de Definition of multivariate distribution families Copyright (C) 2013 -- 2015 Martin Schlather This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suix2te 330, Boston, MA 02111-1307, USA. */ #include "RF.h" #include "Families.h" #include "Operator.h" ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// /* Achtung: falls in *R x gegeben wird, handelt es sich um eine bedingte Verteilung. Da wo NA steht wird simuliert. todo: Fuer *R x koennte man dies auch programmieren, ist aber noch gemacht Falls in *P2sided x[i]==y[i] so wird fuer diese Stellen \frac{\partial}{\prod_{diese i's} \partial y_{i_k}} F(y[j_1], y_{i_k}, ....) - F(x[j_1], y_{i_k}, ....) betrachtet toDo : erweiterung auf asymmetrische Kurven; dann wird nicht mehr in der Mitte aufgeschnitten, sondern wo linke und rechte dichtefunktion sich schneiden */ void arcsqrtP(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double *v) { double scale = 4.0 * P0(ARCSQRT_SCALE), y = *x / scale; if (y <= PIHALF) *v = 0.0; else { *v = atan(sqrt( y / PIHALF - 1.0)) / PIHALF; // without (5) !! // siehe tcf*KS-1.pdf } } // --- not corrected yet: void arcsqrtD(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double *v) { double scale = 4.0 * P0(ARCSQRT_SCALE), y = *x / scale; if (y <= M_PI_2 ) *v = 0.0; else { *v = SQRT2 / (M_PI * scale * y * sqrt(y / M_PI_2 - 2.0)); } } void arcsqrtDlog(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double *v) { arcsqrtD(x, cov, v); *v = log(*v); } void arcsqrtDinverse(double VARIABLE_IS_NOT_USED *v, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *left, double VARIABLE_IS_NOT_USED *right) { if (v == NULL || *v <=0 ) { left[0] = 0; right[0] = RF_INF; } else ERR("Dinverse of arcsqrt unknown"); } void arcsqrtQ(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double *v) { if (*x < 0 || *x > 1) {*v = RF_NA; return;} double scale = 4.0 * P0(ARCSQRT_SCALE); double z = tan(PIHALF * *x); *v = PIHALF * (z * z + 1.0) * scale; } void arcsqrtR(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double *v) { if (x != NULL) *v = *x; else { double u = UNIFORM_RANDOM; arcsqrtQ(&u, cov, v); // printf("scale %f %f \n", scale, *v); assert(*v >= PIHALF); // *v = 0.1; } } int init_arcsqrt(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){ if (cov->mpp.moments >= 0) { cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; } // maxheight based on a Gaussian shape function with the same // sd's but normed to 1 at the origin cov->mpp.maxheights[0] = RF_NA; return NOERROR; } void do_arcsqrt(cov_model *cov, double *v){ arcsqrtR(NULL, cov, v); } void range_arcsqrt(cov_model VARIABLE_IS_NOT_USED *cov, range_type* range){ range->min[ARCSQRT_SCALE] = 0.0; range->max[ARCSQRT_SCALE] = RF_INF; range->pmin[ARCSQRT_SCALE] = 0.0; range->pmax[ARCSQRT_SCALE] = 100000; range->openmin[ARCSQRT_SCALE] = false; range->openmax[ARCSQRT_SCALE] = true; } ////////////////////////////////////////////////////////////////////// void addVariable(char *name, double *x, int nrow, int ncol, SEXP env) { SEXP Y; int j, size= nrow * ncol; if (ncol == 1) PROTECT(Y = allocVector(REALSXP, nrow)); else PROTECT(Y = allocMatrix(REALSXP, nrow, ncol)); //printf("name = %s %f %d %d %d\n",name,x[0],nrow,ncol,isEnvironment(env)); for (j=0; jsexp; int size, nkappas = CovList[cov->nr].kappas, i = nkappas; // PMI(cov); if (cov->ownkappanames != NULL) { while(cov->ownkappanames[--i] != NULL) { addVariable(cov->ownkappanames[i], P(i), cov->nrow[i], cov->ncol[i], env); } } // PMI(cov); // print("which = %d\n", which); // eval( ((sexp_type*) P(1])->sexp , env); res = eval(PLANG(which)->sexp, env); size = P0INT(DISTR_NROW) * P0INT(DISTR_NCOL); for (i=0; isexp); evaluateDistr(cov, DISTR_DX, v); } void distrDlog(double *x, cov_model *cov, double *v) { distrD(x, cov, v); *v = log(*v); } void distrDinverse(double VARIABLE_IS_NOT_USED *v, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED *y) { // v is here INPUT variable NotProgrammedYet("distrDinverse"); } void distrP(double *x, cov_model *cov, double *v) { addVariable((char*) "q", x, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_PX, v); } void distrP2sided(double *x, double *y, cov_model *cov, double *v) { int i, size = P0INT(DISTR_NROW) * P0INT(DISTR_NCOL), dim = cov->xdimown; double w; assert(dim == cov->tsdim); assert(dim == size); if (dim == 1) { double z = x != NULL ? *x : -*y; addVariable((char*) "q", &z, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_PX, &w); addVariable((char*) "q", y, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_PX, v); *v -= w; } else { NotProgrammedYet("multivariate families of distribution functions"); double Sign = 1.0; ALLOC_EXTRA(z, size); for (i=0; i 1) {*v = RF_NA; return;} addVariable((char*) "p", x, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_QX, v); } void distrR(double *x, cov_model *cov, double *v) { if (x != NULL) ERR("Conditional distribution not allowed yet"); addVariable((char*) "n", &ONE, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_RX, v); } void distrR2sided(double *x, double *y, cov_model *cov, double *v) { if (x != NULL || y != NULL) ERR("conditional distribution not allowed yet"); addVariable((char*) "n", &ONE, 1, 1, PENV(DISTR_ENV)->sexp); evaluateDistr(cov, DISTR_RX, v); } void kappa_distr(int i, cov_model *cov, int *nr, int *nc){ // int dim = cov->tsdim; *nc = *nr = i < CovList[cov->nr].kappas ? SIZE_NOT_DETERMINED : -1; } int check_distr(cov_model *cov) { ROLE_ASSERT(ROLE_DISTR); kdefault(cov, DISTR_NROW, 1); kdefault(cov, DISTR_NCOL, 1); cov->vdim[0] = P0INT(DISTR_NROW); cov->vdim[1] = P0INT(DISTR_NCOL); EXTRA_STORAGE; return NOERROR; } int init_distr(cov_model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s){ int err = NOERROR; cov->mpp.maxheights[0] = RF_NA; cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; return err; } void do_distr_do(cov_model *cov, double *v){ distrR(NULL, cov, v); } void range_distr(cov_model *cov, range_type *range){ #define distr_n 5 int idx[distr_n] = {DISTR_DX, DISTR_PX, DISTR_QX, DISTR_RX, DISTR_ENV}; for (int i=0; imin[j] = RF_NAN; range->max[j] = RF_NAN; range->pmin[j] = RF_NAN; range->pmax[j] = RF_NAN; range->openmin[j] = false; range->openmax[j] = false; } range->min[DISTR_NROW] = 1; range->max[DISTR_NROW] = 10; range->pmin[DISTR_NROW] = 1; range->pmax[DISTR_NROW] = 10; range->openmin[DISTR_NROW] = false; range->openmax[DISTR_NROW] = true; range->min[DISTR_NCOL] = 1; range->max[DISTR_NCOL] = 10; range->pmin[DISTR_NCOL] = 1; range->pmax[DISTR_NCOL] = 10; range->openmin[DISTR_NCOL] = false; range->openmax[DISTR_NCOL] = false; int i, kappas = CovList[cov->nr].kappas; for (i=DISTR_LAST + 1; imin[i] = RF_NEGINF; range->max[i] = RF_INF; range->pmin[i] = 1e10; range->pmax[i] = -1e10; range->openmin[i] = true; range->openmax[i] = true; } } #define GAUSS_PARAMETER_BASICS \ double *m = P(GAUSS_DISTR_MEAN), \ *sd = P(GAUSS_DISTR_SD) \ #define GAUSS_PARAMETERS \ GAUSS_PARAMETER_BASICS; \ int i, mi, si,\ len_mean = cov->nrow[GAUSS_DISTR_MEAN], \ len_sd = cov->nrow[GAUSS_DISTR_SD], \ dim = cov->xdimown; \ assert(dim == cov->tsdim); #define FOR for(mi=si=i=0; i 1) {*v = RF_NA; return;} GAUSS_PARAMETER_BASICS; int returnlog = P0INT(GAUSS_DISTR_LOG); *v = qnorm(x[0], m[0], sd[0], true, returnlog); } void gaussR(double *x, cov_model *cov, double *v) { GAUSS_PARAMETERS; if (x == NULL) { FOR v[i] = rnorm(m[mi], sd[si]); } else { FOR v[i] = R_FINITE(x[i]) ? x[i] : rnorm(m[mi], sd[si]); } } void gaussR2sided(double *x, double *y, cov_model *cov, double *v) { GAUSS_PARAMETERS; if (x == NULL) { FOR { while (fabs(v[i] = rnorm(m[mi], sd[si])) > y[i]); } } else { FOR { while ((v[i] = rnorm(m[mi], sd[si])) < x[i] || v[i] > y[i]); } } } void kappa_gauss_distr(int i, cov_model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ if (i == GAUSS_DISTR_SD || i== GAUSS_DISTR_MEAN) { *nc = 1; *nr = SIZE_NOT_DETERMINED; } else *nc = *nr = i == GAUSS_DISTR_LOG ? 1 : -1; } int check_gauss_distr(cov_model *cov) { ROLE_ASSERT(ROLE_DISTR); GAUSS_PARAMETER_BASICS; int dim = cov->xdimown; assert(dim == cov->tsdim); if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; if (m == NULL) kdefault(cov, GAUSS_DISTR_MEAN, 0.0); if (sd == NULL) kdefault(cov, GAUSS_DISTR_SD, 1.0); kdefault(cov, GAUSS_DISTR_LOG, false); cov->vdim[0] = cov->xdimprev; cov->vdim[1] = 1; return NOERROR; } int init_gauss_distr(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){ GAUSS_PARAMETERS; if (cov->mpp.moments >= 0) { cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; if (cov->mpp.moments >= 1) { if (dim > 1) SERR("multivariate moment cannot be calculated"); cov->mpp.mM[1] = cov->mpp.mMplus[1] = m[0]; if (cov->mpp.moments >= 2) { double SD = sd == NULL ? 1.0 : sd[0]; cov->mpp.mM[2] = cov->mpp.mMplus[2] = SD * SD + m[0] * m[0] ; } } } // maxheight based on a Gaussian shape function with the same // sd's but normed to 1 at the origin cov->mpp.maxheights[0] = intpow(SQRTTWOPI, dim); FOR cov->mpp.maxheights[0] *= sd[si]; cov->mpp.unnormedmass = RF_NA; // / cov->mpp.maxheights[0]; cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; return NOERROR; } void do_gauss_distr(cov_model *cov, double *v){ double *sd = P(GAUSS_DISTR_SD); int i, mi, si, len_mean = cov->nrow[GAUSS_DISTR_MEAN], len_sd = cov->nrow[GAUSS_DISTR_SD], dim = cov->xdimown; assert(dim == cov->tsdim); cov->mpp.maxheights[0] = intpow(SQRTTWOPI, -dim); FOR cov->mpp.maxheights[0] /= sd[si]; gaussR(NULL, cov, v); } void range_gauss_distr(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[GAUSS_DISTR_MEAN] = RF_NEGINF; range->max[GAUSS_DISTR_MEAN] = RF_INF; range->pmin[GAUSS_DISTR_MEAN] = -1e8; range->pmax[GAUSS_DISTR_MEAN] = 1e8; range->openmin[GAUSS_DISTR_MEAN] = true; range->openmax[GAUSS_DISTR_MEAN] = true; range->min[GAUSS_DISTR_SD] = 0; range->max[GAUSS_DISTR_SD] = RF_INF; range->pmin[GAUSS_DISTR_SD] = 1e-10; range->pmax[GAUSS_DISTR_SD] = 1e8; range->openmin[GAUSS_DISTR_SD] = true; range->openmax[GAUSS_DISTR_SD] = true; range->min[GAUSS_DISTR_LOG] = 0; range->max[GAUSS_DISTR_LOG] = 1; range->pmin[GAUSS_DISTR_LOG] = 0; range->pmax[GAUSS_DISTR_LOG] = 1; range->openmin[GAUSS_DISTR_LOG] = false; range->openmax[GAUSS_DISTR_LOG] = false; } #define SPHERIC_SPACEDIM 0 #define SPHERIC_BALLDIM 1 #define SPHERIC_RADIUS 2 double random_spheric(int tsdim, int balldim) { int d; double r2 = -1.0; while (r2 < 0.0) { r2 = 1.0; for (d=tsdim; d 1) {*v = RF_NA; return;} ERR("density of 'RRspheric' cannot be calculated yet"); } void sphericR(double *x, cov_model *cov, double *v) { int dim = P0INT(SPHERIC_SPACEDIM), balldim = P0INT(SPHERIC_BALLDIM); if (x==NULL) *v = random_spheric(dim, balldim) * P0(SPHERIC_RADIUS); // q[>0] ist E else ERR("conditional distribution cannot be calculated for sphericP."); } int check_RRspheric(cov_model *cov) { int err; ROLE_ASSERT(ROLE_DISTR); kdefault(cov, SPHERIC_SPACEDIM, 1); kdefault(cov, SPHERIC_BALLDIM, P0INT(SPHERIC_SPACEDIM)); kdefault(cov, SPHERIC_RADIUS, 1.0); if ((err = checkkappas(cov)) != NOERROR) return err; if (cov->tsdim != 1) SERR("only dimension 1 allowed"); cov->vdim[0] = cov->xdimprev; cov->vdim[1] = 1; return NOERROR; } void range_RRspheric(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[SPHERIC_SPACEDIM] = 1; range->max[SPHERIC_SPACEDIM] = RF_INF; range->pmin[SPHERIC_SPACEDIM] = 0; range->pmax[SPHERIC_SPACEDIM] = 10; range->openmin[SPHERIC_SPACEDIM] = false; range->openmax[SPHERIC_SPACEDIM] = true; range->min[SPHERIC_BALLDIM] = 1; range->max[SPHERIC_BALLDIM] = RF_INF; range->pmin[SPHERIC_BALLDIM] = 0; range->pmax[SPHERIC_BALLDIM] = 10; range->openmin[SPHERIC_BALLDIM] = false; range->openmax[SPHERIC_BALLDIM] = true; range->min[SPHERIC_RADIUS] = 0; range->max[SPHERIC_RADIUS] = RF_INF; range->pmin[SPHERIC_RADIUS] = 0.001; range->pmax[SPHERIC_RADIUS] = 1000; range->openmin[SPHERIC_RADIUS] = true; range->openmax[SPHERIC_RADIUS] = true; } int init_RRspheric(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s) { int i, m, nm = cov->mpp.moments, nmvdim = nm + 1, // vdim == 1 dim = P0INT(SPHERIC_SPACEDIM), balldim = P0INT(SPHERIC_BALLDIM), testn = GLOBAL.mpp.n_estim_E; double scale, dummy, *M = cov->mpp.mM, radius = P0(SPHERIC_RADIUS), *Mplus = cov->mpp.mMplus; // printf("ball mppM2 %f\n", cov->mpp.mM2); //printf("%d %d %d\n", testn, GLOBAL.mpp.n_estim_E, // // P0INT(SPHERIC_BALLDIM));// assert(false); for (M[0] = 1.0, m=1; m < nmvdim; M[m++] = 0.0); for (i=0; i 1) PRINTF("init_spheric %f %f %f\n", M[nm], exp( (balldim - dim) * M_LN_SQRT_PI // log(sqrt(pi)) + lgammafn(0.5 * cov->tsdim + 1) - lgammafn(0.5 * balldim + 1)), exp( - dim * M_LN_SQRT_PI // log(sqrt(pi)) + lgammafn(0.5 * cov->tsdim + 1)) ); //cov->mpp.refradius = RF_NA; cov->mpp.maxheights[0] = RF_NA; cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; return NOERROR; } void do_RRspheric(cov_model *cov, double *v) { sphericR(NULL, cov, v); } // ***** deterministic **** #define DETERM_MEAN 0 #define DETERM_PARAMETERS \ double *mean = P(DETERM_MEAN); \ int i, mi, \ dim = cov->xdimown, \ len_mean = cov->nrow[DETERM_MEAN] #define DETERMFOR for(mi=i=0; i mean[mi] || y[i] < mean[mi]) { *v = 0.0; return; } } } else { DETERMFOR { if (x[i] > mean[mi] || y[i] < mean[mi]) { *v = 0.0; return; } } } *v = 1.0; } void determQ(double *x, cov_model *cov, double *v) { if (*x < 0 || *x > 1) {*v = RF_NA; return;} v[0] = P0(DETERM_MEAN); } void determR(double *x, cov_model *cov, double *v) { DETERM_PARAMETERS; if (x==NULL) DETERMFOR v[i] = mean[i]; else DETERMFOR v[i] = !R_FINITE(x[i]) || x[i] == mean[mi] ? mean[mi] : RF_NA; } void kappa_determ(int i, cov_model *cov, int *nr, int *nc){ int dim = cov->tsdim; *nc = 1; *nr = i == 0 ? dim : i == 1 ? 1 : -1; } void determR2sided(double *x, double *y, cov_model *cov, double *v) { DETERM_PARAMETERS; if (x == NULL) { DETERMFOR v[i] = fabs(y[i]) > mean[mi] ? mean[mi] : RF_NA; } else { DETERMFOR v[i] = x[i] < mean[mi] && y[i] > mean[mi] ? mean[mi] : RF_NA; } } int check_determ(cov_model *cov) { int dim = cov->xdimown; if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; if (PisNULL(DETERM_MEAN)) kdefault(cov, DETERM_MEAN, 0.0); cov->vdim[0] = dim; cov->vdim[1] = 1; return NOERROR; } void range_determ(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[DETERM_MEAN] = RF_NEGINF; range->max[DETERM_MEAN] = RF_INF; range->pmin[DETERM_MEAN] = -1e10; range->pmax[DETERM_MEAN] = 1e10; range->openmin[DETERM_MEAN] = true; range->openmax[DETERM_MEAN] = true; } int init_determ(cov_model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) { int err = NOERROR; // normed against an extremal gaussian process with shape that // is identically constant 1 in space cov->mpp.maxheights[0] = 1.0; cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; return err; } void do_determ(cov_model *cov, double *v) { determR(NULL, cov, v); } // ***** Set **** void setParamD(double *x, cov_model *cov, double *v) { VTLG_D(x, cov->sub[SETPARAM_LOCAL], v); } void setParamDlog(double *x, cov_model *cov, double *v) { VTLG_DLOG(x, cov->sub[SETPARAM_LOCAL], v); } void setParamDinverse(double *v, cov_model *cov, double *left, double *right) { NONSTATINVERSE_D(v, cov->sub[SETPARAM_LOCAL], left, right); } void setParamP(double *x, cov_model *cov, double *v) { VTLG_P(x, cov->sub[SETPARAM_LOCAL], v); } void setParamP2sided(double *x, double *y, cov_model *cov, double *v) { VTLG_P2SIDED(x, y, cov->sub[SETPARAM_LOCAL], v); } void setParamQ(double *x, cov_model *cov, double *v) { VTLG_Q(x, cov->sub[SETPARAM_LOCAL], v); } void setParamR(double *x, cov_model *cov, double *v) { VTLG_R(x, cov->sub[SETPARAM_LOCAL], v); } void setParamR2sided(double *x, double *y, cov_model *cov, double *v) { VTLG_R2SIDED(x, y, cov->sub[SETPARAM_LOCAL], v); } int check_setParam(cov_model *cov) { cov_model *next= cov->sub[SETPARAM_LOCAL]; int err, dim = cov->xdimown; kdefault(cov, SET_PERFORMDO, true); if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; if ((err = CHECK_R(next, dim)) != NOERROR) return err; setbackward(cov, next); cov->vdim[0] = next->vdim[0]; cov->vdim[1] = next->vdim[1]; TaylorCopy(cov, next); cov->mpp.maxheights[0] = next->mpp.maxheights[0]; cov->mpp.unnormedmass = next->mpp.unnormedmass; return NOERROR; } int init_setParam(cov_model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) { cov_model *next= cov->sub[SETPARAM_LOCAL]; set_storage *X = cov->Sset; // PMI(cov->calling); crash(); assert(X != NULL); int err; if ((err = INIT(next, cov->mpp.moments, s)) != NOERROR) return err; if (X->remote != NULL) { assert(X->set != NULL); X->set(cov->sub[0], X->remote, X->variant); } TaylorCopy(cov, next); cov->mpp.maxheights[0] = next->mpp.maxheights[0]; cov->mpp.unnormedmass = next->mpp.unnormedmass; return NOERROR; } void do_setParam(cov_model *cov, double *v) { cov_model *next = cov->sub[SETPARAM_LOCAL]; bool performDo = P0INT(SET_PERFORMDO); if (performDo) { DORANDOM(next, v); } //cov->mpp.maxheights[0] = next->mpp.maxheights[0]; } void range_setParam(cov_model *cov, range_type *range){ range_setparam(cov, range); /* range->min[SETPARAM_VARIANT] = 0; range->max[SETPARAM_VARIANT] = 0; range->pmin[SETPARAM_VARIANT] = 0; range->pmax[SETPARAM_VARIANT] = 0; range->openmin[SETPARAM_VARIANT] = false; range->openmax[SETPARAM_VARIANT] = false; */ } #define LOC_PARAMETER_BASICS \ cov_model *next = cov->sub[0]; \ int dim=cov->xdimown; \ double *loc = P(LOC_LOC), \ *scale= P(LOC_SCALE) #define LOC_PARAMETERS \ LOC_PARAMETER_BASICS; \ int i, mi, si, \ len_loc = cov->nrow[LOC_LOC], \ len_scale = cov->nrow[LOC_SCALE]; \ assert(dim == cov->tsdim) #define LOCFOR for(mi=si=i=0; i 1.0 / 100.0 && *v != 1.0) { // LOCFOR { // printf("loc %d %f %f s=%f %f\n", i, y[i], z[i], scale[si], loc[mi]); // } // printf("v=%e\n", *v); //} } else { ALLOC_DOLLAR2(zy, dim); LOCFOR { z[i] = (x[i] - loc[mi]) / scale[si]; zy[i] = (y[i] - loc[mi]) / scale[si]; } VTLG_P2SIDED(z, zy, next, v); } } void locQ(double *x, cov_model *cov, double *v) { LOC_PARAMETER_BASICS; if (dim != 1) BUG; VTLG_Q(x, next, v); v[0] = v[0] * scale[0] + loc[0]; } void locR(double *x, cov_model *cov, double *v) { LOC_PARAMETERS; double *z1 = NULL; if (x != NULL) { ALLOC_DOLLAR(Z1, dim); z1 = Z1; LOCFOR z1[i] = (x[i] - loc[mi]) / scale[si]; } VTLG_R(z1, next, v); if (x == NULL) LOCFOR v[i] = v[i] * scale[si] + loc[mi]; else LOCFOR v[i] = R_FINITE(x[i]) ? x[i] : v[i] * scale[si] + loc[mi]; // printf("v=%f %f %f s=%f l=%f\n", v[0], v[1], v[2], scale[0], loc[0]); // APMI(cov->calling); } void locR2sided(double *x, double *y, cov_model *cov, double *v) { LOC_PARAMETERS; double *z1 = NULL; if (x != NULL) { ALLOC_DOLLAR(Z1, dim); z1 = Z1; LOCFOR z1[i] = (x[i] - loc[mi]) / scale[si]; } // !! Z1 can be != NULL as used differently among the loc-fcts ALLOC_DOLLAR2(z2, dim); assert(y != NULL); LOCFOR z2[i] = (y[i] - loc[mi]) / scale[si]; VTLG_R2SIDED(z1, z2, next, v); // printf("locR dim=%d scale=%f %f %f %f -> %f %f %f %d -> %f %f %f\n", dim, *scale, y[0], y[1], y[2], z2[0], z2[1], z2[2], z1==NULL, v[0], v[1], v[2]); LOCFOR { v[i] = v[i] * scale[si] + loc[mi]; } } void kappa_loc(int i, cov_model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ if (i == LOC_SCALE || i== LOC_LOC) { *nc = 1; *nr = SIZE_NOT_DETERMINED; } else if (i == LOC_POWER) *nc = * nr = 1; else *nc = *nr = -1; } int check_loc(cov_model *cov) { ROLE_ASSERT(ROLE_DISTR); cov_model *next = cov->sub[0]; int // len_loc = cov->nrow[LOC_LOC], // len_scale = cov->nrow[LOC_SCALE], dim = cov->xdimown; if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; double *loc = P(LOC_LOC), *scale = P(LOC_SCALE); int err; kdefault(cov, POWPOWER, 0.0); //PMI(cov); if ((err = CHECK_R(next, dim)) != NOERROR) return err; // assert(false); /// not checked yet if (loc == NULL) kdefault(cov, LOC_LOC, 0.0); if (scale == NULL) kdefault(cov, LOC_SCALE, 1.0); cov->vdim[0] = next->vdim[0]; cov->vdim[1] = next->vdim[1]; DOLLAR_STORAGE; return NOERROR; } int init_loc(cov_model *cov, gen_storage *s){ // cov_fct *C = CovList + cov->nr; LOC_PARAMETERS; int err; double p = P0(LOC_POWER); // same role as POWPOWER of '$power' if ((err = INIT(next, cov->mpp.moments, s)) != NOERROR) return err; if (cov->mpp.moments >= 0) { cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; if (cov->mpp.moments >= 1) { if (dim > 1) { LOCFOR if (scale[si] != 1.0 || loc[mi]!=0.0) SERR("multivariate moment cannot be calculated"); } cov->mpp.mM[1] = cov->mpp.mM[1] * scale[0] + loc[0]; cov->mpp.mMplus[1] = loc[0] == 0.0 ? cov->mpp.mMplus[1] * scale[0] : RF_NA; if (cov->mpp.moments >= 2) { double ssq = scale[0] * scale[0]; cov->mpp.mM[2] = cov->mpp.mM[2] * ssq + loc[0] * (2.0 * cov->mpp.mM[1] - loc[0]); cov->mpp.mMplus[1] = loc[0] == 0.0 ? cov->mpp.mMplus[1] * ssq : RF_NA; } } } // normed against a shape function with the same scale modification if (R_FINITE(next->mpp.unnormedmass)) cov->mpp.unnormedmass = next->mpp.unnormedmass * pow(scale[0], dim + p); else cov->mpp.maxheights[0] = next->mpp.maxheights[0] / scale[0]; cov->mpp.mM[0] = next->mpp.mM[0]; cov->mpp.mMplus[0] = next->mpp.mMplus[0]; return NOERROR; } void do_loc(cov_model *cov, double *v){ // cov_model *next = cov->sub[0]; //double *scale= P(LOC_SCALE); DORANDOM(cov->sub[0], v); locR(NULL, cov, v); // cov->mpp.maxheights[0] = next->mpp.maxheights[0] * scale[0]; } void range_loc(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[LOC_LOC] = RF_NEGINF; range->max[LOC_LOC] = RF_INF; range->pmin[LOC_LOC] = -1e8; range->pmax[LOC_LOC] = 1e8; range->openmin[LOC_LOC] = true; range->openmax[LOC_LOC] = true; range->min[LOC_SCALE] = 0; range->max[LOC_SCALE] = RF_INF; range->pmin[LOC_SCALE] = 1e-10; range->pmax[LOC_SCALE] = 1e8; range->openmin[LOC_SCALE] = true; range->openmax[LOC_SCALE] = true; range->min[LOC_POWER] = RF_NEGINF; range->max[LOC_POWER] = RF_INF; range->pmin[LOC_POWER] = -cov->xdimown; range->pmax[LOC_POWER] = +cov->xdimown; range->openmin[LOC_POWER] = true; range->openmax[LOC_POWER] = true; } #define UNIF_PARAMETER_BASICS \ double \ *min = (double*) P(UNIF_MIN), \ *max= (double*) P(UNIF_MAX) #define UNIF_PARAMETERS \ UNIF_PARAMETER_BASICS; \ int i, \ mini, \ maxi, \ len_min = cov->nrow[UNIF_MIN], \ len_max = cov->nrow[UNIF_MAX], \ dim = cov->xdimown; \ assert(dim == cov->tsdim) #define UNIFOR for(mini=maxi=i=0; i max[maxi]) { *v = 0.0; return; } else if (normed) area *= max[maxi] - min[mini]; } *v = 1.0 / area; } void unifDlog(double *x, cov_model *cov, double *v) { unifD(x, cov, v); *v = log(*v); } void unifDinverse(double *v, cov_model *cov, double *left, double *right) { UNIF_PARAMETERS; double area = 1.0; if (P0INT(UNIF_NORMED)) UNIFOR area *= max[maxi] - min[mini]; if (*v * area > 1){ UNIFOR left[i] = right[i] = 0.5 * (max[maxi] + min[mini]); } else { UNIFOR { left[i] = min[mini]; right[i] = max[maxi]; // printf("unif=%d %f %f\n", i, min[mini], max[maxi]); } } } void unifP(double *x, cov_model *cov, double *v) { // distribution functions UNIF_PARAMETERS; double area = 1.0; bool normed = P0INT(UNIF_NORMED); UNIFOR { if (x[i] <= min[mini]) { *v = 0.0; return; } if (x[i] < max[maxi]) area *= x[i] - min[mini]; if (normed) area /= max[maxi] - min[mini]; // ekse area *= 1.0; } *v = area; } void unifP2sided(double *x, double *y, cov_model *cov, double *v) { // distribution functions UNIF_PARAMETERS; double a, b; bool normed = P0INT(UNIF_NORMED); double area = 1.0; UNIFOR { a = x != NULL ? x[i] : -y[i]; if (a != y[i]) { if (a < min[mini]) a = min[mini]; b = y[i] > max[maxi] ? max[maxi] : y[i]; if (a >= b) { *v = 0.0; return; } area *= b - a; } else { if (a < min[mini] || a > max[maxi]) { *v = 0.0; return; } } if (normed) area /= max[maxi] - min[mini]; } *v = area; } void unifQ(double *x, cov_model *cov, double *v) { if (*x < 0 || *x > 1) {*v = RF_NA; return;} UNIF_PARAMETER_BASICS; if (P0INT(UNIF_NORMED)) { *v = min[0] + *x * (max[0] - min[0]); } else { *v = min[0] + *x; } } void unifR(double *x, cov_model *cov, double *v) { UNIF_PARAMETERS; if (x == NULL) { UNIFOR { v[i] = min[mini] + UNIFORM_RANDOM * (max[maxi] - min[mini]); } } else { UNIFOR { v[i] = !R_FINITE(x[i]) ? min[mini] + UNIFORM_RANDOM * (max[maxi] - min[mini]) : x[i] >= min[mini] && x[i] <= max[maxi] ? x[i] : RF_NA; } } } void unifR2sided(double *x, double *y, cov_model *cov, double *v) { UNIF_PARAMETERS; double a, b; UNIFOR { a = x != NULL ? (x[i] < min[mini] ? min[mini] : x[i]) : (-y[i] < min[mini] ? min[mini] : -y[i]); b = y[i] > max[maxi] ? max[maxi] : y[i]; if (a > b) ERR("parameters of 2-sided unifR out of range"); v[i] = a + UNIFORM_RANDOM * (b-a); } } void kappa_unif(int i, cov_model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ if (i == UNIF_MIN || i == UNIF_MAX) { *nc = 1; *nr = SIZE_NOT_DETERMINED; } else if (i == UNIF_NORMED) { *nc = *nr = 1; } else *nc = *nr = -1; } int check_unif(cov_model *cov) { //PMI(cov->calling); ROLE_ASSERT(ROLE_DISTR); int // len_min = cov->nrow[UNIF_MIN], // len_max = cov->nrow[UNIF_MAX], dim = cov->xdimown; if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; if (PisNULL(UNIF_MIN)) kdefault(cov, UNIF_MIN, 0.0); if (PisNULL(UNIF_MAX)) kdefault(cov, UNIF_MAX, 1.0); kdefault(cov, UNIF_NORMED, true); cov->vdim[0] = cov->tsdim; cov->vdim[1] = 1; //PMI(cov->calling); return NOERROR; } int init_unif(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){ UNIF_PARAMETERS; // printf("entering init_unif\n"); // normed against a shape function that is the indicator function // of the respective window cov->mpp.unnormedmass = 1.0; UNIFOR { cov->mpp.unnormedmass *= max[maxi] - min[mini]; // printf("mini=%d %d; %f %f %f\n", mini, maxi, max[maxi], min[mini], // cov->mpp.unnormedmass); } if (P0INT(UNIF_NORMED)) { cov->mpp.maxheights[0] = 1.0 / cov->mpp.unnormedmass; if (cov->mpp.moments >= 0) { cov->mpp.mM[0] = cov->mpp.mMplus[0] = 1.0; if (cov->mpp.moments >= 1) { if (dim > 1) SERR("multivariate moment cannot be calculated"); cov->mpp.mM[1] = 0.5 * (min[0] + max[0]); cov->mpp.mMplus[1] = max[0] <= 0.0 ? 0.0 : 0.5 * max[0]; if (cov->mpp.moments >= 2) { cov->mpp.mM[2] = max[0] - min[0]; cov->mpp.mM[2] *= cov->mpp.mM[2] / 12.0; } } } } else { // not normed cov->mpp.maxheights[0] = 1.0; cov->mpp.mM[0] = cov->mpp.mMplus[0] = cov->mpp.unnormedmass; if (cov->mpp.moments > 0) SERR("unnormed unif does not allow for higher moments"); } // PMI(cov); // printf("leaving init_unif\n"); return NOERROR; } void do_unif(cov_model *cov, double *v){ unifR(NULL, cov, v); // printf("v=%f\n", *v); // cov->mpp.maxheights[0] = 1.0; // UNIFOR cov->mpp.maxheights[0] /= max[maxi] - min[mini]; } void range_unif(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[UNIF_MIN] = RF_NEGINF; range->max[UNIF_MIN] = RF_INF; range->pmin[UNIF_MIN] = -1e8; range->pmax[UNIF_MIN] = 1e8; range->openmin[UNIF_MIN] = true; range->openmax[UNIF_MIN] = true; range->min[UNIF_MAX] = RF_NEGINF; range->max[UNIF_MAX] = RF_INF; range->pmin[UNIF_MAX] = -1e8; range->pmax[UNIF_MAX] = 1e8; range->openmin[UNIF_MAX] = true; range->openmax[UNIF_MAX] = true; range->min[UNIF_NORMED] = 0; range->max[UNIF_NORMED] = 1; range->pmin[UNIF_NORMED] = 0; range->pmax[UNIF_NORMED] = 1; range->openmin[UNIF_NORMED] = false; range->openmax[UNIF_NORMED] = false; } #define SURFACE(dim) ((dim) * intpow(2.0, dim)) //of a d-dim cube with edge length 2, i.e. a (dim-1) dimensional surface ! // Ring constanter intensitaet #define TwoDSimu \ X = UNIFORM_RANDOM * (start + end); \ Y = (2.0 * UNIFORM_RANDOM - 1.0) * (end - start); \ i0 = UNIFORM_RANDOM < 0.5; \ \ x[1 - i0] = Y >= 0 ? Y + start : Y - start; \ x[i0] = ((Y >= 0) xor (i0==1)) ? X - start : start - X; double VolumeOfCubeRing(double *xsort, double start, double end, int dim, int squeezed_parts) { int d, red_dim = dim - squeezed_parts; assert(red_dim > 0); double res = intpow(2.0, dim); // printf("res = %f\n", res); for (d=1; d<=squeezed_parts; d++) res *= xsort[d]; // printf("res = %f %f %d %f\n", res, xsort[1], red_dim, // intpow(end, red_dim) - intpow(start, red_dim)); return res * (intpow(end, red_dim) - intpow(start, red_dim)); } double PoweredVolOfCube(double *xsort, double start, double end, double p, int dim, int squeezed_parts) { double red_dim = dim - squeezed_parts; assert(red_dim > 0); double res = red_dim * intpow(2.0, dim), // surface * 2^squeezed pPd = p + red_dim; int d; for (d=1; d<=squeezed_parts; d++) res *= xsort[d]; //printf("res=%f %f %f %f\n", res, pow(end, pPd), pow(start, pPd), pPd); //return * intpow(2 * start, squeezed_parts) return res * (pow(end, pPd) - pow(start, pPd)) / pPd; } double ExpVolOfCube(double start, double end, double p, double s, int dim, int squeezed_parts) { assert(start >= 0 && end > start); // aus evaluate_rectangular // die Dichte gleich c * s * p * x^{p-1-(d-1)} * e^{-s x^p} / b_d // Fuer nicht squeezed Dimensionen d' = d - squeezed_parts gilt (mit c=1) // int_start^end c * s * p * x^{p-d} * e^{-s x^p} / b_d d x // = // b_{d'} * int_start^end c * s * p * x^{p-d} * x^{d'-1} * e^{-s x^p}/b_d d x // =_{u = s x^p; du= s p x^{p-1} dx} // b_{d'}/b_d * int_{s start^p}^{s end^p} * x^{d'-d} * e^{-u} d u // = // b_{d'}/b_d * s^{(d - d')/p} * // int_{s start^p}^{s end^p} * u^{(d'-d) / p} * e^{-u} d u double red_dim = dim - squeezed_parts, a = - squeezed_parts / p + 1.0; assert(red_dim > 0); // printf("s=%f start=%f p=%f sq=%d a=%f %f\n", s, start, p, squeezed_parts, a, s * pow(start, p)); // printf("expvol %f %f %f a=%f ret=%f\n", // SURFACE(red_dim) / SURFACE(dim), pow(s, squeezed_parts / p), // incomplete_gamma(s * pow(start, p), s * pow(end, p), a), a, // SURFACE(red_dim) / SURFACE(dim) * pow(s, squeezed_parts / p) * // incomplete_gamma(s * pow(start, p), s * pow(end, p), a)); return SURFACE(red_dim) / SURFACE(dim) * pow(s, squeezed_parts / p) * incomplete_gamma(s * pow(start, p), s * pow(end, p), a); } void RandomPointOnCubeRing(double start, double end, int dim, double *x) { double twostart = 2.0 * start, X, Y; int i0; // printf("start = %f %f\n", start, end); assert(start >= 0 && end > start); // todo linear on the ring instead of constant ?! // todo general dimension switch(dim) { case 1: { x[0] = (2.0 * UNIFORM_RANDOM - 1.0) * (end - start); if (x[0] < 0.0) x[0] -= start; else x[0] += start; } break; case 2: TwoDSimu break; case 3: { // Aufsplitten in Roehre mit obigen 2-d Schnitt und Laenge 2*start // und den beiden noch fehlenden Kappen der Groesse 4 end^2 double massTube, massEnds, Z; massTube = 4.0 * (start + end) * (end - start) * twostart; massEnds = 2.0 * end; massEnds = 2.0 * massEnds * massEnds; assert(massTube > 0); assert(massEnds > 0); X = UNIFORM_RANDOM * (massTube + massEnds); if (X < massTube) { // Roehre TwoDSimu x[2] = (2.0 * UNIFORM_RANDOM - 1.0) * start; } else { // kappen x[0] = (2.0 * UNIFORM_RANDOM - 1.0) * end; x[1] = (2.0 * UNIFORM_RANDOM - 1.0) * end; Z = (2.0 * UNIFORM_RANDOM - 1.0) * (end - start); x[2] = Z > 0 ? Z + start : Z - start; } //printf("x= %f %f %f start=%f %f Z=%f\n", x[0], x[1],x[2], start, end, Z); } break; default : BUG; } //printf("x0= %f %f %f dim=%d\n", x[0], x[1],x[2],dim); } #define TwoDSurface \ X = UNIFORM_RANDOM * 2 * dist; \ if (X > 4 * dist) { \ if (X > 6 * dist) { \ x[0] = -dist; \ x[1] = X - 7 * dist; \ } else { \ x[1] = dist; \ x[0] = X - 5 * dist; \ } \ } else { \ if (X > 2 * dist) { \ x[0] = dist; \ x[1] = X - 3 * dist; \ } else { \ x[1] = -dist; \ x[0] = X - 1 * dist; \ } \ } void RandomPointOnCubeSurface(double dist, int dim, double *x) { double X; assert(dist >= 0); // todo general dimension // printf("RandomPointOnCubeSurface dim=%d\n", dim); switch(dim) { case 1: x[0] = UNIFORM_RANDOM < 0.5 ? dist : -dist; break; case 2: TwoDSurface break; case 3: { if ( (X = UNIFORM_RANDOM * 6) > 2 ) { // Roehre TwoDSurface x[2] = (2.0 * UNIFORM_RANDOM - 1.0) * dist; } else { // Kappen x[0] = (2.0 * UNIFORM_RANDOM - 1.0) * dist; x[1] = (2.0 * UNIFORM_RANDOM - 1.0) * dist; x[2] = X > 1 ? -dist : dist; } } break; default : BUG; } } // rectangular distribution // Marco's Idee um Kirstins Sachen zu implementieren #define IDX_INNER 0 #define IDX_STEPS 1 #define IDX_OUTER (NSTEP + IDX_STEPS) // change also KeyInfo if changed #define ASSIGN_IDX_INNER (-IDX_STEPS) #define ASSIGN_IDX_OUTER (-IDX_STEPS - 1) #define ASSIGN_IDX_MISMATCH (-IDX_STEPS - 2) #define INNER rect->inner #define INNER rect->inner #define INNER_CONST rect->inner_const #define INNER_POW rect->inner_pow #define OUTER rect->outer #define OUTER_CONST rect->outer_const #define OUTER_POW rect->outer_pow #define OUTER_POW_CONST rect->outer_pow_const #define STEP rect->step #define NSTEP rect->nstep #define VALUE(IDX) rect->value[IDX_STEPS + IDX] #define INNER_CUM rect->weight[IDX_INNER] #define OUTER_CUM rect->weight[IDX_OUTER] #define WEIGHT(IDX) rect->weight[IDX_STEPS + IDX] #define TMP (rect->tmp_n) #define TMP_WEIGHT rect->tmp_weight #define RIGHT_END(IDX) rect->right_endpoint[IDX] #define SQUEEZED(IDX) rect->squeezed_dim[IDX] #define ASSIGNMENT(IDX) rect->asSign[IDX] // for control only #define RECTANGULAR_PARAMETER_BASICS \ rect_storage *rect = cov->Srect; \ int VARIABLE_IS_NOT_USED dim = cov->xdimown; \ assert(dim == cov->tsdim); \ if (rect == NULL) BUG #define RECTANGULAR_PARAMETERS \ RECTANGULAR_PARAMETER_BASICS; \ cov_model VARIABLE_IS_NOT_USED *next = cov->sub[0]; \ #define RECTFOR for(i=0; iweight // bereits bekannt ist ? RECTANGULAR_PARAMETERS; double min, *ysort = rect->ysort, *zneu = rect->z; int d, kstep, dM1, start_cumul = 1, dimP1 = dim + 1, One = 1; bool use_weights = cumsum != rect->weight; // printf("y=%f \n", y[0]); ysort[0] = 0.0; for (d=0; di; if (dim > 1) { RU_ordering(ysort, dimP1, One, i); assert(dim < 3 || (ysort[i[0]] <= ysort[i[1]] && ysort[i[1]] <= ysort[i[2]])); } else { i[0] = 0; i[1] = 1; } assert(i[0] == 0 && i[1] > 0); //for (d=0; d<=dim; d++) printf("%d,%f ", i[d], z[d]); printf("\n"); for(d=0; d= INNER && use_weights) { kstep = zneu[d] >= OUTER ? NSTEP : (int) ((zneu[d] - INNER) / STEP); // if (DEBUG_CUMSUM) printf("kstep %d %f\n", kstep, (zneu[d] - INNER)/STEP);// if (kurz) { if (zneu[d] == RF_INF) { cumsum[TMP] = OUTER_CUM; TMP++; assert(cumsum[TMP-1] >= 0); return; } RIGHT_END(TMP) = INNER + kstep * STEP; SQUEEZED(TMP) = - MAXINT / 2; ASSIGNMENT(TMP)= ASSIGN_IDX_MISMATCH; // if (DEBUG_CUMSUM) printf("0. TMP=%d kstep=%d ass=%d\n", TMP, kstep, ASSIGNMENT(TMP));// cumsum[TMP++] = WEIGHT(kstep - 1); assert(cumsum[TMP-1] >= 0); } else { int k; for (k=-IDX_STEPS; k= 0); } } zneu[dM1] = INNER + STEP * kstep; start_cumul = TMP; } else { kstep = 0; // printf("AA d=%d %d %f %f\n", d, TMP, zneu[dM1], INNER); } for ( ; d<=dim; d++) { dM1 = d - 1; if (zneu[d] == zneu[dM1]) continue; //printf("%f == %f %d\n",zneu[d], zneu[dM1], zneu[d] == zneu[dM1] ); // printf("A d=%d %d\n", d, TMP); if (zneu[dM1] < INNER) { min = RIGHT_END(TMP) = zneu[d] <= INNER ? zneu[d] : INNER; SQUEEZED(TMP) = dM1; ASSIGNMENT(TMP) = ASSIGN_IDX_INNER; cumsum[TMP++] = INNER_CONST * PoweredVolOfCube(ysort, zneu[dM1], min, INNER_POW, dim, dM1); // if (DEBUG_CUMSUM) printf("2. TMP=%d kstep=%d ass=%d cum=%f inner_c=%f start=%f ende=%f pow=%f dim=%d squ=%d\n", TMP, kstep, ASSIGNMENT(TMP), cumsum[TMP-1], // INNER_CONST, zneu[dM1], min,PoweredVolOfCube(ysort, zneu[dM1], min, INNER_POW, dim, dM1), dim, dM1);// assert(false); assert(cumsum[TMP-1] >= 0); if (zneu[d] <= INNER) continue; zneu[dM1] = INNER; } if (zneu[dM1] < OUTER) { min = zneu[d] <= OUTER ? zneu[d] : OUTER; int k, steps = (min - zneu[dM1]) / STEP; double a = zneu[dM1]; assert(steps <= NSTEP); //printf("d=%d min=%f %f %f\n", d, min, zneu[dM1], zneu[d]); // assert(TMP <= NSTEP + 2 + dim); for (k=0; ksub[0], &z); // cumsum[TMP++] = VolumeOfCubeRing(ysort, a, a + STEP, dim, dM1) * VALUE(kstep); assert(cumsum[TMP-1] >= 0); // if (DEBUG_CUMSUM) printf("3. TMP=%d kstep=%d ass=%d a=%4.2f st=%4.2f end=%4.2f val=%4.2e, tot=%4.4f\n", TMP, kstep, ASSIGNMENT(TMP), a, STEP, RIGHT_END(TMP), // VALUE(kstep), cumsum[TMP-1]); // assert(dM1 != 1); } // // printf("B %d %d\n", d, TMP); assert(TMP <= NSTEP + 2 + dim); if (kstep < NSTEP) { RIGHT_END(TMP) = min; SQUEEZED(TMP) = dM1; ASSIGNMENT(TMP) = kstep; assert(TMP > 0); // if (DEBUG_CUMSUM) printf("4. TMP=%d k=%d ass=%d\n", TMP, kstep, ASSIGNMENT(TMP));// cumsum[TMP++] = VolumeOfCubeRing(ysort, a, min, dim, dM1) * VALUE(kstep); assert(cumsum[TMP-1] >= 0); if (zneu[d] <= OUTER) continue; } // printf("outer=%f %f %f %e %d %d %d\n", OUTER, a, cumsum[TMP-1], // OUTER - a, OUTER == min, dM1, kstep < NSTEP); //assert(false); // printf("C %d %d\n", d, TMP); } RIGHT_END(TMP) = zneu[d]; SQUEEZED(TMP) = dM1; ASSIGNMENT(TMP) = ASSIGN_IDX_OUTER; // if (DEBUG_CUMSUM) printf("5. TMP=%d k=%d ass=%d %f p=%f right=%f\n", // TMP, 9999, ASSIGNMENT(TMP), OUTER, OUTER_POW, zneu[d]); if (next->finiterange == true) { cumsum[TMP++] = 0.0; continue; } else if (OUTER_POW > 0) { // printf("here %d C=%f Z=%f PC=%f e=%f %f\n", /// IDX_OUTER, OUTER_CONST,zneu[d], // OUTER_POW_CONST, // ExpVolOfCube(OUTER, zneu[d], OUTER_POW, // OUTER_POW_CONST, dim, dM1), // OUTER_POW_CONST * ExpVolOfCube(OUTER, zneu[d], OUTER_POW, // OUTER_POW_CONST, dim, dM1) // ); cumsum[TMP++] = OUTER_CONST * ExpVolOfCube(OUTER, zneu[d], OUTER_POW, OUTER_POW_CONST, dim, dM1); // if (DEBUG_CUMSUM) printf("6. TMP=%d k=%d ass=%d oc=%f c=%f op=%f opc=%f dim=%d %d right=%f\n", TMP, 10001, ASSIGNMENT(TMP),OUTER_CONST, OUTER, OUTER_POW, // OUTER_POW_CONST, dim, dM1, zneu[d]); // printf("TMP=%d %f\n", TMP, cumsum[TMP-1]); // if (cumsum[TMP-1] < 0) { //t i; //fori=1; i= 0); assert(R_FINITE(cumsum[TMP-1])); } else { //printf("here XXX\n"); cumsum[TMP++]= OUTER_CONST * PoweredVolOfCube(ysort, OUTER, zneu[d], OUTER_POW, dim, dM1); assert(cumsum[TMP-1] >= 0); } } //PMI(cov); // make cummulative weights // //d=0; printf("cumsum %d %e; %d %e\n", d, cumsum[d], start_cumul-1, cumsum[start_cumul-1]); for (d=start_cumul; d= ASSIGNMENT(d-1)); } //APMI(cov); } // (unnormierte) VERTEILUNGSFUNKTION F(x,x,...,x) ist // c * (exp(-O^p - exp(-a x^p)) fuer grosse x, also integriert ueber // den Raum, gesehen als Funktion einer Variable x; O=OUTER // Somit ist die Dichte gleich c * s * p * x^{p-d} * e^{-s x^p} / b // wobei b die Oberflaeche des Wuerfels mit Kantenlaenge 2 bezeichnet // // weiter in ExpVolOfCube #define OUTER_DENSITY(v, x) \ if (OUTER_POW > 0) { \ double pow_x = OUTER_POW_CONST * pow(x, OUTER_POW); \ v = OUTER_CONST * OUTER_POW * pow_x * intpow(x, -dim) * \ exp(-pow_x) / SURFACE(dim); \ } else { \ /* DICHTE ist c*x^p fuer grosse x */ \ v = OUTER_CONST * pow(x, OUTER_POW); \ } void evaluate_rectangular(double *x, cov_model *cov, double *v) { // *x wegen searchInverse !! RECTANGULAR_PARAMETERS; if (*x < 0.0) BUG; //crash(); assert(*x >= 0); if (*x <= INNER) { // for technical reasons <= not <; see Dinverse // DICHTE ist c*x^p nahe Ursprung // printf("inner "); *v = INNER_CONST * pow(*x, INNER_POW); // printf("x=%f %f %f\n", *x, INNER, *v); return; } else if (*x >= OUTER) { //printf("outer "); if (next->finiterange == true) { *v = 0.0; return; } OUTER_DENSITY(*v, *x); return; } else { // printf("k=%d %f %f ", (int) ((*x - INNER) / STEP), // VALUE(0), // VALUE( (int) ((*x - INNER) / STEP))); *v = VALUE( (int) ((*x - INNER) / STEP)); return; } } void rectangularD(double *x, cov_model *cov, double *v) { bool onesided = P0INT(RECT_ONESIDED); if (onesided && *x <= 0) { *v = 0; return; } if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); RECTANGULAR_PARAMETER_BASICS; int i; double max = RF_NEGINF; RECTFOR {double y=fabs(x[i]); if (y > max) max = y;} evaluate_rectangular(&max, cov, v); // printf("max %f %f\n", max, *v); // printf("rectD %f %f %d %f\n", max, *v, P0INT(RECT_NORMED), OUTER_CUM); if (P0INT(RECT_NORMED)) *v /= OUTER_CUM; if (onesided) *v *= 2.0; // R^1 !! } void rectangularDlog(double *x, cov_model *cov, double *v) { rectangularD(x, cov, v); *v = log(*v); } void rectangularDinverse(double *V, cov_model *cov, double *left, double *right) { if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); RECTANGULAR_PARAMETERS; double er, outer, x = RF_NA, v = *V; int i; bool onesided = P0INT(RECT_ONESIDED); if (P0INT(RECT_NORMED)) v *= OUTER_CUM; // unnormed v if (onesided) v /= 2.0; // printf("V=%e dim=%d\n", *V, dim); // assert(*V >= 0.0 && *V <= 1.0); if (*V <= 0.0) { RECTFOR { left[i] = RF_NEGINF; right[i] = RF_INF; //printf("i=%d %f %f\n", i, left[i], right[i]); } return; } if (!next->finiterange == true && OUTER_POW > 1) { // local maximum of density function exists for outer_pow > 1 outer = pow( (OUTER_POW -1) / (OUTER_POW * OUTER_POW_CONST), 1 / OUTER_POW); if (outer < OUTER) outer = OUTER; } else outer = OUTER; evaluate_rectangular(&outer, cov, &er); // printf("outer =%f %f %f %f\n", outer, OUTER, er, v); if (er > v) { // inverse is in the tail // first guess: double inverse, eps = 0.01; if (OUTER_POW > 0) { // printf("here\n"); // rough evaluation inverse = pow(- log(v / (OUTER_CONST * OUTER_POW)) / OUTER_POW_CONST, 1.0 / OUTER_POW); if (inverse <= outer) inverse = 2 * outer; x = searchInverse(evaluate_rectangular, cov, inverse, outer, v, eps); } else // printf("hereB\n"); x = pow(OUTER_CONST / v, 1 / OUTER_POW); } else { //printf("hereC\n"); int nsteps = (OUTER - INNER) / STEP; for (i=nsteps - 1; i>=0; i--) if (VALUE(i) >= v) { x = INNER + (i+1) * STEP; break; } if (i<0) { // all values had been smaller, so the // inverse is in the inner part // printf("hereD\n"); evaluate_rectangular(&(INNER), cov, &er); if (er >= v) x = INNER; else if (INNER_POW == 0) x = 0.0; else if (INNER_POW < 0.0) { x = pow(v / INNER_CONST, 1.0 / INNER_POW); } else BUG; } } RECTFOR { //printf("rect i=%d x=%f\n", x); left[i] = onesided ? 0.0 : -x; right[i] = x; } } void rectangularP(double VARIABLE_IS_NOT_USED *x, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *v) { if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); // RECTANGULAR_PARAMETERS; NotProgrammedYet(""); } void rectangularP2sided(double *x, double *y, cov_model *cov, double *v) { // distribution functions bool onesided = P0INT(RECT_ONESIDED); int d; if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); RECTANGULAR_PARAMETER_BASICS; if (x != NULL) BUG; if (onesided && *y <= 0) { *v = 0; return; } for (d=0; d= 0.0); if (y[d] == 0.0) { *v = 0.0; return; } } if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); // CumSum(y, false, cov, rect->tmp_weight); // double tw = TMP_WEIGHT[TMP-1]; CumSum(y, true, cov, rect->tmp_weight); *v = TMP_WEIGHT[TMP-1]; // print("tw=%f %f y=%f %f %f\n", tw, *v, y[0], y[1], y[2]); // if(*v != tw) APMI(cov); // assert(*v == tw); if (P0INT(RECT_NORMED)) *v /= OUTER_CUM; // printf("p2sided %f %d outercum=%f, %d %f\n", // *y, P0INT(RECT_NORMED), OUTER_CUM, // TMP, TMP_WEIGHT[TMP-1]); // int i; for (i=0; i 1) {*v = RF_NA; return;} if (!P0INT(RECT_APPROX)) ERR("approx=FALSE only for simulation"); NotProgrammedYet("rectangularQ"); } #define ACCEPT_REJECT \ if (!P0INT(RECT_APPROX)) { \ int ni = dim, \ quot = dim + 1, \ bytes = sizeof(double) * dim; \ double newquot, approx, truevalue, \ max = RF_NEGINF; \ RECTFOR {double z=fabs(v[i]); if (z > max) max = z;} \ evaluate_rectangular(&max, cov, &approx); \ ABSFCTN(v, next, &truevalue); \ newquot = truevalue / approx; \ assert(quot < cov->qlen + 2); \ if (isMonotone(next->monotone)) { /* rejection sampling */ \ assert(approx >= truevalue); \ cov->q[ni] = 0.0; \ if (UNIFORM_RANDOM >= newquot) { \ /*printf(".//");*/ \ RESTART; \ } /* else printf("://"); */ \ } else { /* MCMC */ \ if (R_FINITE(cov->q[ni])) { \ cov->q[ni] -= - 1.0; \ if (UNIFORM_RANDOM * cov->q[quot] >= newquot) { /* keep old one */ \ memcpy(v, cov->q, bytes); \ } else { /* accept new one*/ \ cov->q[quot] = newquot; /* pi / Q */ \ memcpy(cov->q, v, bytes); \ } \ } else { /* !R_FINITE(cov->q[ni]), i.e. starting value */ \ cov->q[ni] = P0INT(RECT_MCMC_N) - 1.0; \ cov->q[quot] = newquot; /* pi / Q */ \ memcpy(cov->q, v, bytes); \ } \ } \ if (cov->q[ni] > 0) { \ RESTART; \ } else { \ cov->q[ni] = P0INT(RECT_MCMC_N); \ /* if (final) { */ \ /* cov->total_n++;*/ \ /* cov->total_sum += truevalue;*/ \ /*} */ \ return; \ } \ } else { \ if (final) { \ double approx, \ max = RF_NEGINF; \ RECTFOR {double z=fabs(v[i]); if (z > max) max = z;} \ evaluate_rectangular(&max, cov, &approx); \ /* cov->total_n++; */ \ /* cov->total_sum += approx; */ \ } \ } void rectangularR(double *x, cov_model *cov, double *v) { //printf("rectR %lu\n", x); if (x != NULL) { //crash(); ERR("put 'flat = false'"); } RECTANGULAR_PARAMETERS; bool final = true; EntryPoint_R: int i = CeilIndex(UNIFORM_RANDOM * OUTER_CUM, rect->weight, NSTEP + 2); //int i = searchFirstGreater(rect->weight, NSTEP + 2, UNIFORM_RANDOM * OUTER_CUM); //printf("i=%d dim=%d %d %d %s\n", i, dim, IDX_INNER, IDX_OUTER, NICK(cov->sub[0])); if (i == IDX_INNER) { RandomPointOnCubeSurface(pow(UNIFORM_RANDOM, 1.0 / (INNER_POW+dim)) * INNER, dim, v); } else if (i == IDX_OUTER) { double u; assert(next->finiterange == false); // if (OUTER_POW > 0) // sampling by inversion formula u = pow(pow(OUTER, OUTER_POW) - log(UNIFORM_RANDOM) / OUTER_POW_CONST, 1 / OUTER_POW); else u = pow(UNIFORM_RANDOM, 1.0 / (OUTER_POW + dim)) * OUTER; // p+d < 0 ! RandomPointOnCubeSurface(u, dim, v); } else { // i >= 0 i -= IDX_STEPS; double start = INNER + i * STEP; RandomPointOnCubeRing(start, start + STEP, dim, v); } if (P0INT(RECT_ONESIDED)) *v = fabs(*v); #define RESTART goto EntryPoint_R ACCEPT_REJECT; } static int zi=0, zs=0, zo=0; void rectangularR2sided(double *x, double *y, cov_model *cov, double *v) { if (x != NULL) NotProgrammedYet("2-sided distribution function for rectangular"); RECTANGULAR_PARAMETERS; int d, sel, red_dim, i, *I = rect->i; double *u, start, end, *ysort = rect->ysort; EntrancePoint_R2: CumSum(y, false, cov, rect->tmp_weight); // ACHTUNG! reck->z gesetzt //assert(rect->tmp_weight[TMP-1] >= 2.1 && rect->tmp_weight[TMP-1] <= 2.2); double random = UNIFORM_RANDOM * rect->tmp_weight[TMP-1]; bool final = SQUEEZED(TMP-1) == 0 && (!P0INT(RECT_APPROX) || !next->deterministic); // ohne Einschraenkung assert(!final); // nur zur Kontrolle sel = CeilIndex(random, rect->tmp_weight, TMP); //printf("r=%f inner=%f o.tmp=%f lst.step.tmp=%f\n", random, rect->tmp_weight[IDX_INNER], rect->tmp_weight[TMP-1],rect->tmp_weight[TMP-1]); red_dim = dim - SQUEEZED(sel); if (red_dim <= 0) BUG; start = sel > 0 ? RIGHT_END(sel - 1) : 0; end = RIGHT_END(sel); assert(start < end); // printf("%f dim=%d squ=%d sel=%d, %f %f %d\n", // y[0], dim, SQUEEZED(sel), sel, random, rect->tmp_weight[TMP-1], TMP); u = rect->tmp_weight; // nutzen des Platzes... //printf("assess %d %d %d\n", ASSIGNMENT(sel), sel, TMP); if (ASSIGNMENT(sel) == ASSIGN_IDX_INNER) { // RandomPointOnCubeSurface(pow(UNIFORM_RANDOM, 1.0 / (INNER_POW+dim)) // * INNER, dim, v); zi ++; double // a + b x^p probab distr on [s, e], i.e. b = 1/(e^p-s^p), p = INNER_POW + red_dim, sp = pow(start, p), ep = pow(end, p), binv = ep - sp, a = - sp / binv, r = pow((UNIFORM_RANDOM - a) * binv, 1 / p); // // printf("inner %d %d ap=%f r=%f inner=%f w=%f tmp=%f, %f\n", dim, SQUEEZED(sel), ap, r, INNER, INNER_CUM, TMP_WEIGHT[IDX_INNER], TMP_WEIGHT[IDX_INNER+1]); RandomPointOnCubeSurface(r, red_dim, u); } else if (ASSIGNMENT(sel) == ASSIGN_IDX_OUTER) { zo++; // printf("outer\n"); double w; assert(next->finiterange == false);// if (OUTER_POW > 0) {// sampling by inversion formula double op = pow(OUTER, OUTER_POW), factor = 1.0 - exp(- OUTER_POW_CONST * (pow(end, OUTER_POW) - op)); w = pow(op - log(1.0 - UNIFORM_RANDOM * factor) / OUTER_POW_CONST, 1 / OUTER_POW); } else w =pow(1.0 - UNIFORM_RANDOM *(1.0 - pow(end/OUTER, OUTER_POW + red_dim)), 1.0 / (OUTER_POW + red_dim)) * OUTER; // p+d < 0 ! //printf("hXXere\n"); RandomPointOnCubeSurface(w, red_dim, u); } else { // i \in steps zs++; // printf("steps \n"); assert(ASSIGNMENT(sel) >= 0 && ASSIGNMENT(sel) < NSTEP); // printf("hAere\n"); // APMI(cov->calling); RandomPointOnCubeRing(start, end, red_dim, u); } // printf("zi=%d s=%d o=%d tot=%d\n", zi, zs, zo, zi+zs+zo); // up to now simulation on cubes on the non-squeezed dimensions // Now the coordinates have to put on the right place and // the remaining coordinates have to be filled with uniformly scattered // random variables on [-ysort[d], ysort[d] // I[.] and ysort[.] have been set by CumSum assert(I[0] == 0.0); for (d=1; d<=SQUEEZED(sel); d++) { v[I[d] - 1] = (2.0 * UNIFORM_RANDOM - 1.0) * ysort[d];//I[d] indiziert ab 1 // da ysort[0]=0 per def. // printf("%f\n", ysort[d]); // v[I[d] - 1] = 2 * UNIFORM_RANDOM - 1; /* print */ } for (d=SQUEEZED(sel); dtaylor[0][TaylorPow]; bool inpownot0 = in_pow != 0.0; assert(!PisNULL(RECT_MAXSTEPS)); int i, j, max_steps = P0INT(RECT_MAXSTEPS), parts = P0INT(RECT_PARTS), max_iterations = P0INT(RECT_MAXIT); assert(next->taylor[0][TaylorConst] > 0.0); INNER_CONST = next->taylor[0][TaylorConst] * safetyP1; INNER_POW = inpownot0 ? in_pow * EMsafety - dimsafety : 0.0; // INNER i = 0; INNER = starting_point; ABSFCTN(&(INNER), next, &v); m = INNER_CONST * pow(INNER, INNER_POW); // printf("inner=%e m=%e v=%e inner_c=%4.3f (%4.3f) inner_p=%4.3f (%4.3f)\n", INNER, m, v, INNER_CONST, next->taylor[0][TaylorConst], INNER_POW, in_pow); while (m < v && INNER >= inner_min) { INNER *= inner_factor; INNER_CONST *= safetyP1; if (inpownot0) INNER_POW = INNER_POW * EMsafety - dimsafety; ABSFCTN(&(INNER), next, &v); m = INNER_CONST * pow(INNER, INNER_POW); // printf("inner=%e m=%e v=%e i=%d, maxit=%d\n", // INNER, m, v, i, max_iterations); } if (INNER < inner_min && m < v) { //PMI(cov); SERR("no majorant found close to the origin"); } i = 0; while (true) { old_ratio = m / v; delta = INNER / parts; double x = INNER - delta; //printf("delta %f x=%f i=%d\n", delta, x, i); for (j=1; j= old_ratio || (!inpownot0 && v <= INNER_CONST) )) old_ratio = ratio; } if (j == parts) break; INNER *= inner_factor; if (INNER > x) INNER = x; if (++i > max_iterations) SERR2("%d iterations performed without success. Increase the value of '%s'", max_iterations, distr[RECT_MAXIT]); } //printf("inner j=%d (%d) inner=%4.3f m=%4.3f\n", // j, max_iterations, INNER, m); // OUTER if (next->tail[0][TaylorExpPow] <= 0.0) { // polynomial tail OUTER_POW = next->tail[0][TaylorPow] * safetyP1 + dimsafety; OUTER_POW_CONST = RF_NA; } else { // exponential tail assert(next->tail[0][TaylorExpConst] > 0.0); OUTER_POW = next->tail[0][TaylorExpPow] / safetyP1; OUTER_POW_CONST = next->tail[0][TaylorExpConst] / safetyP1; } OUTER_CONST = next->tail[0][TaylorConst] * safetyP1; if (next->finiterange == true) { // printf("xx %s %d\n", CovList[next->sub[0]->nr].nick, CovList[next->sub[0]->nr].finiterange == 1); // APMI(cov); INVERSE(ZERO, next, &(OUTER)); if (INNER >= OUTER) OUTER = INNER * (1.0 + 1e-5); } else { // ! finiterange assert(next->tail[0][TaylorConst] > 0.0); i = 0; OUTER = starting_point * safetyP1; ABSFCTN(&(OUTER), next, &v); // m = OUTER_CONST * (OUTER_POW <= 0 ? pow(OUTER, OUTER_POW) // : exp(- OUTER_POW_CONST * pow(OUTER, OUTER_POW))); OUTER_DENSITY(m, OUTER); while (m < v && OUTER < outer_max) { if (++i > max_iterations) SERR1("No majorant found. Function does not allow for a majorant or increase '%s'", distr[RECT_MAXIT]); OUTER_CONST *= safetyP1; OUTER_POW /= safetyP1; OUTER *= outer_factor; ABSFCTN(&(OUTER), next, &v); OUTER_DENSITY(m, OUTER); // printf("i =%d outer %f m=%e v=%e\n", i, OUTER, m, v); assert(false); } if (OUTER > outer_max && m < v) SERR("No majorant found close for large arguments"); // printf("i =%d outer %f m=%e v=%e\n", i, OUTER, m, v); //assert(false); //printf("i =%d outer %f %f m=%f v=%f\n", //i, OUTER, // OUTER_CONST * (OUTER_POW <= 0 ? pow(OUTER, OUTER_POW) // : exp(- OUTER_POW_CONST * pow(OUTER, OUTER_POW))), // m , v); //assert(false); i=0; while (true) { old_ratio = m / v; delta = OUTER / parts; double x = OUTER + delta; for (j=1; j max_iterations) BUG; } } // printf("inner outer %f %f\n", INNER, OUTER);// assert(false); // STEPS // todo variable Schrittweiten waeren wegen des steilen // Anstiegs am Pol besser STEP = (OUTER - INNER) / max_steps; if (STEP < min_steps) { NSTEP = (int) ((OUTER - INNER) / min_steps); STEP = (OUTER - INNER) / NSTEP; } else NSTEP = max_steps; // PMI(cov, "getm"); return NOERROR; } int check_rectangular(cov_model *cov) { cov_model *next = cov->sub[0]; int err, dim = cov->xdimown; ASSERT_CARTESIAN; ROLE_ASSERT(ROLE_DISTR); kdefault(cov, RECT_SAFETY, GLOBAL.distr.safety); kdefault(cov, RECT_MINSTEPLENGTH, GLOBAL.distr.minsteplen); kdefault(cov, RECT_MAXSTEPS, GLOBAL.distr.maxsteps); kdefault(cov, RECT_PARTS, GLOBAL.distr.parts); kdefault(cov, RECT_MAXIT, GLOBAL.distr.maxit); kdefault(cov, RECT_INNERMIN, GLOBAL.distr.innermin); kdefault(cov, RECT_OUTERMAX, GLOBAL.distr.outermax); kdefault(cov, RECT_MCMC_N, GLOBAL.distr.mcmc_n); kdefault(cov, RECT_NORMED, true); kdefault(cov, RECT_APPROX, true); kdefault(cov, RECT_ONESIDED, false); // int onesided = P0INT(RECT_ONESIDED); if (cov->q == NULL) QALLOC(dim + 2); cov->q[dim] = RF_NA; if ((err = CHECK(next, dim, dim, ShapeType, XONLY, dim == 1 && P0INT(RECT_ONESIDED) ? CARTESIAN_COORD : ISOTROPIC, SCALAR, ROLE_DISTR)) != NOERROR) { return err; } if (!next->deterministic) { //APMI(cov); SERR("currently, only deterministic submodels are allowed"); // to do } // if (!isMonotone(next->monotone)) SERR("only monotone submodels are allowed"); if (next->taylorN <= 0 || next->tailN <= 0) { //PMI(next); SERR1("'%s' does not have integrability information", NICK(next)); } assert(R_FINITE(next->taylor[0][TaylorPow])); if (next->taylor[0][TaylorPow] <= -dim) SERR1("pole of '%s' not integrable", NICK(next)); assert(R_FINITE(next->tail[0][TaylorPow])); // PMI(next); if (next->tail[0][TaylorPow] >= -dim && next->tail[0][TaylorExpPow] == 0.0 && next->tail[0][TaylorConst] != 0.0) SERR1("tail of '%s' not integrable", NICK(next)); assert(R_FINITE(next->tail[0][TaylorConst])); //PMI(next); if (next->taylor[0][TaylorConst] == 0.0) { //PMI(next); SERR1("'%s' seems to be a trivial shape function", NICK(next)); } // APMI(cov); if (cov->xdimprev != dim || dim != cov->tsdim) return ERRORDIM; cov->vdim[0] = cov->tsdim; cov->vdim[1] = 1; //PMI(cov); return NOERROR; } int init_rectangular(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){ int err; NEW_STORAGE(rect); RECTANGULAR_PARAMETERS; // muss nach obigen stehen und vor allem anderen if ((err = INIT(next, cov->mpp.moments, s)) != NOERROR) { //PMI(cov); return err; } if ((err = GetMajorant(cov)) != NOERROR) return err; // muss genau nach MALLOC oben und MALLOC unten stehen if (INNER >= OUTER) BUG; // APMI(cov); int d, abschnitte = NSTEP + 2, tmp_n = NSTEP + 2 + dim; double x; if ((rect->value = (double*) MALLOC(sizeof(double) * abschnitte)) == NULL || (rect->weight = (double*) MALLOC(sizeof(double) * abschnitte)) == NULL|| (rect->tmp_weight = (double*) CALLOC(tmp_n, sizeof(double))) == NULL || (rect->right_endpoint = (double*) MALLOC(sizeof(double) * tmp_n)) ==NULL|| (rect->ysort = (double*) MALLOC(sizeof(double) * (dim + 1))) == NULL || (rect->z = (double*) MALLOC(sizeof(double) * (dim + 1))) == NULL ||//dummy (rect->squeezed_dim = (int*) MALLOC(sizeof(int) * tmp_n)) == NULL || (rect->asSign = (int*) MALLOC(sizeof(int) * tmp_n)) == NULL || (rect->i = (int*) MALLOC(sizeof(int) * (dim + 1))) == NULL) return ERRORMEMORYALLOCATION; x = INNER; for (d=0; d= 0); } rect->value[IDX_INNER] = rect->value[IDX_OUTER] = RF_NA; for(d=0; dtmp_weight[d] = RF_INF; CumSum(rect->tmp_weight, false, cov, rect->weight); cov->mpp.mM[0] = cov->mpp.mMplus[0] = P0INT(RECT_NORMED) ? 1.0 : OUTER_CUM; assert(R_FINITE(OUTER_CUM)); if (cov->mpp.moments > 0) { cov->mpp.mM[1] = next->mpp.mM[1]; cov->mpp.mMplus[1] = next->mpp.mMplus[1]; if (!R_FINITE(cov->mpp.mM[1])){ // crash(); APMI(cov); BUG; } } // normed against cov->sub[0] cov->mpp.maxheights[0] = RF_NA; // cov->mpp.unnormedmass = OUTER_CUM; // if (OUTER_CUM < 1.0) { APMI(cov); } if (false) { double mini=1e-6, u,w,t, tot_u=0.0, tot_t=0.0, e, xx = mini; // mini = STEP; xx=INNER-mini; double end; int k, n = NSTEP ; n = 1 * NSTEP; assert(mini >= 0.0); assert(xx >= 0.0); for (k=0, end=INNER; k<=n+1; k++, end+=STEP) { if (k==n+1) end *= 5; u = t = 0.0; while (xx < end) { ABSFCTN(&xx, next, &w); double y = xx + 1e-10; assert(y >= 0.0); evaluate_rectangular(&y, cov, &e); //printf("%f %e %e\n", xx, w, e); u += w * 4 * M_PI / 3 * (powl(xx + mini, 3) - powl(xx, 3)); t += w * (powl(2 * (xx + mini), 3) - powl(2 * xx, 3)); xx += mini; } tot_u += u; tot_t += t; // printf("%d %e %e tot=%f %f\n", k-1, u, t, tot_u, tot_t); } } // printf("totalweight=%f\n", OUTER_CUM); assert(OUTER_CUM >= 1.0 || next->nr != STROKORB_MONO); // // printf("init rect %d %d\n", NSTEP, TMP); assert(NSTEP == TMP - 2); //cov->total_n = 0; // cov->total_sum = 0.0; return NOERROR; } void do_rectangular(cov_model *cov, double *v){ cov_model *next = cov->sub[0]; //int err; gen_storage s; gen_NULL(&s); DO(next, &s); rectangularR(NULL, cov, v); } void range_rectangular(cov_model *cov, range_type *range){ range->min[RECT_SAFETY] = 0; range->max[RECT_SAFETY] = RF_INF; range->pmin[RECT_SAFETY] = 0.001; range->pmax[RECT_SAFETY] = 0.1; range->openmin[RECT_SAFETY] = true; range->openmax[RECT_SAFETY] = true; range->min[RECT_MINSTEPLENGTH] = 0; range->max[RECT_MINSTEPLENGTH] = RF_INF; range->pmin[RECT_MINSTEPLENGTH] = 0.01; range->pmax[RECT_MINSTEPLENGTH] = 10; range->openmin[RECT_MINSTEPLENGTH] = false; range->openmax[RECT_MINSTEPLENGTH] = true; range->min[RECT_MAXSTEPS] = 1; range->max[RECT_MAXSTEPS] = RF_INF; range->pmin[RECT_MAXSTEPS] = 2; range->pmax[RECT_MAXSTEPS] = 10000; range->openmin[RECT_MAXSTEPS] = false; range->openmax[RECT_MAXSTEPS] = true; range->min[RECT_PARTS] = 2; range->max[RECT_PARTS] = RF_INF; range->pmin[RECT_PARTS] = 3; range->pmax[RECT_PARTS] = 10; range->openmin[RECT_PARTS] = false; range->openmax[RECT_PARTS] = true; range->min[RECT_MAXIT] = 1; range->max[RECT_MAXIT] = RF_INF; range->pmin[RECT_MAXIT] = 2; range->pmax[RECT_MAXIT] = 10; range->openmin[RECT_MAXIT] = false; range->openmax[RECT_MAXIT] = true; range->min[RECT_INNERMIN] = 0; range->max[RECT_INNERMIN] = RF_INF; range->pmin[RECT_INNERMIN] = 1e-100; range->pmax[RECT_INNERMIN] = 1; range->openmin[RECT_INNERMIN] = true; range->openmax[RECT_INNERMIN] = true; range->min[RECT_OUTERMAX] = 0; range->max[RECT_OUTERMAX] = RF_INF; range->pmin[RECT_OUTERMAX] = 1; range->pmax[RECT_OUTERMAX] = 1e10; range->openmin[RECT_OUTERMAX] = true; range->openmax[RECT_OUTERMAX] = true; range->min[RECT_MCMC_N] = 1; range->max[RECT_MCMC_N] = MAXINT; range->pmin[RECT_MCMC_N] = 0; range->pmax[RECT_MCMC_N] = 1000; range->openmin[RECT_MCMC_N] = false; range->openmax[RECT_MCMC_N] = true; range->min[RECT_NORMED] = 0; range->max[RECT_NORMED] = 1; range->pmin[RECT_NORMED] = 0; range->pmax[RECT_NORMED] = 1; range->openmin[RECT_NORMED] = false; range->openmax[RECT_NORMED] = false; range->min[RECT_APPROX] = 0; range->max[RECT_APPROX] = 1; range->pmin[RECT_APPROX] = 0; range->pmax[RECT_APPROX] = 1; range->openmin[RECT_APPROX] = false; range->openmax[RECT_APPROX] = false; range->pmin[RECT_ONESIDED] = range->min[RECT_ONESIDED] = 0; range->pmax[RECT_ONESIDED] = range->max[RECT_ONESIDED] = cov->tsdim == 1; range->openmin[RECT_ONESIDED] = false; range->openmax[RECT_ONESIDED] = false; } #define MCMC_FCTN 0 #define MCMC_MCMC_N 0 #define MCMC_MCMC_SIGMA 1 #define MCMC_NORMED 2 #define MCMC_MAXDENSITY 3 #define MCMC_RANDLOC 4 #define MCMC_GIBBS 5 #define MCMC_LAST_PARAM MCMC_GIBBS void kappa_mcmc(int i, cov_model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ *nc = 1; *nr = (i == MCMC_MCMC_SIGMA) ? 0 : (i <= MCMC_LAST_PARAM) ? 1 : -1; } void mcmcIntegral(cov_model VARIABLE_IS_NOT_USED *cov) { NotProgrammedYet("mcmcIntegral"); /* int n = ; bool normed = P0INT(MCMC_NORMED); double value; PINT(MCMC_NORMED)[0] = false; while (cov->q[MCMC_CURR_N] >= 1.0) { cov->q[MCMC_CURR_SUM] -= 1.0; mcmcR(NULL, unklar welche Funktion, &value); } */ } void mcmcD(double *x, cov_model *cov, double *v) { cov_model *next = cov->sub[MCMC_FCTN]; // int mcmc_n = P0INT(MCMC_MCMC_N); FCTN(x, next, v); *v = fabs(*v); if (P0INT(MCMC_NORMED)) { BUG; // not programmed yet } } void mcmcDlog(double *x, cov_model *cov, double *v) { mcmcD(x, cov, v); *v = log(*v); } void mcmcDinverse(double VARIABLE_IS_NOT_USED *V, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *left, double VARIABLE_IS_NOT_USED *right) { NotProgrammedYet("mcmcDinverse"); } void mcmcP(double VARIABLE_IS_NOT_USED *x, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *v) { NotProgrammedYet("mcmcP"); } void mcmcP2sided(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED *y, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *v) { NotProgrammedYet("mcmcP2sided"); } void mcmcQ(double *x, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *v) { if (*x < 0 || *x > 1) {*v = RF_NA; return;} NotProgrammedYet("mcmcQ"); } void mcmcR(double *x, cov_model *cov, double *v) { if (x != NULL) { ERR("put 'flat = false'"); } cov_model *next = cov->sub[MCMC_FCTN]; location_type *loc=Loc(cov); int i, d, vdim = cov->tsdim, n = P0INT(MCMC_MCMC_N); double proposedvalue, maxdens = P0(MCMC_MAXDENSITY), *sigma = P(MCMC_MCMC_SIGMA), posvalue = cov->Smcmc->posvalue, *delta = cov->Smcmc->deltapos, *pos = cov->Smcmc->pos; assert(delta != NULL && pos != NULL); bool gibbs = (bool) P0INT(MCMC_GIBBS), randloc = (bool) P0INT(MCMC_RANDLOC); ALLOC_NEW(Smcmc, proposed, vdim, proposed); ALLOC_NEW(Smcmc, propdelta, vdim, propdelta); for (i=0; inrow[MCMC_MCMC_SIGMA]])); } else { for (d=0; d < vdim; d++) { // printf("%d %d %f %f\n", d, d % cov->nrow[MCMC_MCMC_SIGMA], // sigma[d % cov->nrow[MCMC_MCMC_SIGMA]], proposed[d]); proposed[d] = (propdelta[d] += rnorm(0.0, sigma[d % cov->nrow[MCMC_MCMC_SIGMA]])); } } if (loc != NULL && randloc) { // to do if (loc->grid) { for (d = 0; d < vdim; d++) { proposed[d] += loc->xgr[d][XSTART] + loc->xgr[d][XSTEP] * (double) ((int) UNIFORM_RANDOM * (loc->xgr[d][XLENGTH] - 1.0)); } } else { double *xx = loc->x + vdim * (int) (loc->spatialtotalpoints*UNIFORM_RANDOM); if (loc->Time) { int vM1 = vdim - 1; for (d = 0; d < vM1; d++) proposed[d] += xx[d]; proposed[d] += loc->T[XSTART] + loc->T[XSTEP] * (double) ((int) UNIFORM_RANDOM * (loc->T[XLENGTH] - 1.0)); } else { for (d = 0; d < vdim; d++) proposed[d] += xx[d]; } } } FCTN(proposed, next, &proposedvalue); if (proposedvalue > maxdens) { proposedvalue = maxdens; } if (proposedvalue > posvalue || proposedvalue > posvalue*UNIFORM_RANDOM) { posvalue = proposedvalue; for (d = 0; d < vdim; d++) { pos[d] = proposed[d]; delta[d] = propdelta[d]; } } } // for iSmcmc->posvalue = posvalue; for (d = 0; d < vdim; d++) v[d] = pos[d]; } void mcmcR2sided(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED *y, cov_model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *v) { } int check_mcmc(cov_model *cov) { cov_model *next = cov->sub[MCMC_FCTN]; int err, vdim; ASSERT_CARTESIAN; ROLE_ASSERT(ROLE_DISTR); kdefault(cov, MCMC_NORMED, false); if (P0INT(MCMC_NORMED)) NotProgrammedYet("mcmc (normed=TRUE)"); vdim = cov->tsdim; if (vdim != cov->xdimown) SERR("inconsistent dimensions given."); if ((err = CHECK(next, vdim, vdim, ShapeType, XONLY, CARTESIAN_COORD, SCALAR, ROLE_DISTR)) != NOERROR) { return err; } cov->vdim[0] = vdim; cov->vdim[1] = 1; if (PisNULL(MCMC_MCMC_SIGMA)) { location_type *loc = Loc(next); if (loc != NULL && loc->grid) { PALLOC(MCMC_MCMC_SIGMA, vdim, 1); for (int d=0; dxgr[d][XSTEP] * 0.1; } else { SERR1("'%s' must be given.", KNAME(MCMC_MCMC_SIGMA)); } } kdefault(cov, MCMC_MCMC_N, GLOBAL.distr.mcmc_n); kdefault(cov, MCMC_MAXDENSITY, 1000); kdefault(cov, MCMC_RANDLOC, false); // scattering among // the locations if on a grid? kdefault(cov, MCMC_GIBBS, false); NEW_STORAGE(mcmc); return NOERROR; } int init_mcmc(cov_model *cov, gen_storage VARIABLE_IS_NOT_USED *s){ cov_model *next = cov->sub[MCMC_FCTN]; location_type *loc=Loc(cov); int d, err, vdim = cov->tsdim; double maxdens = P0(MCMC_MAXDENSITY); if ((err = INIT(next, cov->mpp.moments, s)) != NOERROR) return err; ALLOC_NEW(Smcmc, pos, vdim, pos); ALLOC_NEW(Smcmc, delta, vdim, deltapos); //PMI(cov); for (d=0; dlx > 0) { if (loc->grid) { assert(loc->xgr[0]!= NULL); for (d=0; dxgr[d][XSTART]; } else if (loc->Time) { int vM1 = vdim - 1; assert(loc->x != NULL); for (d=0; dx[d]; pos[d] = loc->T[XSTART]; } else { assert(loc->x != NULL); for (d=0; dx[d]; } } FCTN(pos, next, &(cov->Smcmc->posvalue)); if (cov->Smcmc->posvalue > maxdens) cov->Smcmc->posvalue = maxdens; return NOERROR; } void do_mcmc(cov_model *cov, double *v){ cov_model *next = cov->sub[MCMC_FCTN];//int err; gen_storage s; gen_NULL(&s); DO(next, &s); mcmcR(NULL, cov, v); } void range_mcmc(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[MCMC_MCMC_N] = 1; range->max[MCMC_MCMC_N] = RF_INF; range->pmin[MCMC_MCMC_N] = 20; range->pmax[MCMC_MCMC_N] = 50; range->openmin[MCMC_MCMC_N] = false; range->openmax[MCMC_MCMC_N] = true; range->min[MCMC_MCMC_SIGMA] = 0; range->max[MCMC_MCMC_SIGMA] = RF_INF; range->pmin[MCMC_MCMC_SIGMA] = 1e-7; range->pmax[MCMC_MCMC_SIGMA] = 1e7; range->openmin[MCMC_MCMC_SIGMA] = true; range->openmax[MCMC_MCMC_SIGMA] = true; range->min[MCMC_NORMED] = 0; range->max[MCMC_NORMED] = 1; range->pmin[MCMC_NORMED] = 0; range->pmax[MCMC_NORMED] = 1; range->openmin[MCMC_NORMED] = false; range->openmax[MCMC_NORMED] = false; range->min[MCMC_MAXDENSITY] = 0; range->max[MCMC_MAXDENSITY] = RF_INF; range->pmin[MCMC_MAXDENSITY] = 1e-7; range->pmax[MCMC_MAXDENSITY] = 1e7; range->openmin[MCMC_MAXDENSITY] = true; range->openmax[MCMC_MAXDENSITY] = true; range->min[MCMC_RANDLOC] = 0; range->max[MCMC_RANDLOC] = 1; range->pmin[MCMC_RANDLOC] = 0; range->pmax[MCMC_RANDLOC] = 1; range->openmin[MCMC_RANDLOC] = false; range->openmax[MCMC_RANDLOC] = false; range->min[MCMC_GIBBS] = 0; range->max[MCMC_GIBBS] = 1; range->pmin[MCMC_GIBBS] = 0; range->pmax[MCMC_GIBBS] = 1; range->openmin[MCMC_GIBBS] = false; range->openmax[MCMC_GIBBS] = false; }