21#ifndef JKQTPSTATKDE_H_INCLUDED
22#define JKQTPSTATKDE_H_INCLUDED
37#include "jkqtmath/jkqtmath_imexport.h"
38#include "jkqtmath/jkqtplinalgtools.h"
39#include "jkqtmath/jkqtparraytools.h"
40#include "jkqtcommon/jkqtpdebuggingtools.h"
41#include "jkqtmath/jkqtpstatbasics.h"
70 return exp(-0.5*fabs(t))/2.0;
78 return (fabs(t)<1.0)?(0.75*(1.0-t*t)):0.0;
86 return (fabs(t)<=1.0)?0.5:0.0;
94 return (fabs(t)<=1.0)?(1.0-fabs(t)):0.0;
103 return (fabs(t)<=1.0)?(15.0/16.0*
jkqtp_sqr(1.0-t*t)):0.0;
111 return (fabs(t)<1.0)?(35.0/32.0*
jkqtp_cube(1.0-t*t)):0.0;
159 return (fabs(tx)<1.0 && fabs(ty)<=1.0)?0.25:0.0;
191template <
class InputIt>
195 return 1.06*sigma/pow(
static_cast<double>(N), 1.0/5.0);
214template <
class InputIt>
215inline double jkqtpstatEvaluateKernelSum(
double t, InputIt first, InputIt last,
const std::function<
double(
double)>& kernel,
double bandwidth) {
218 for (
auto it=first; it!=last; ++it) {
221 const double vx=(t-v)/bandwidth;
226 if (cnt==0)
return 0.0;
227 return res/
static_cast<double>(cnt)/bandwidth;
254template <
class InputIt,
class OutputIt>
255inline void jkqtpstatKDE1DAutoranged(InputIt first, InputIt last, OutputIt KDEXOut, OutputIt KDEYOut,
int Nout=100,
const std::function<
double(
double)>& kernel=std::function<
double(
double)>(&
jkqtpstatKernel1DGaussian),
double bandwidth=1.0,
bool cummulative=
false) {
256 double minV=0, maxV=0;
258 jkqtpstatMinMax<InputIt>(first, last, minV, maxV,
nullptr,
nullptr, &N);
260 std::vector<double> histX;
261 std::vector<double> histY;
263 const double range=maxV-minV;
264 const double binw=range/
static_cast<double>(Nout);
267 for (
double xi=minV; xi<=maxV; xi+=binw) {
271 if (histX.size()>0 && histX[histX.size()-1]<maxV) {
272 histX.push_back(maxV);
280 for (
size_t i=0; i<histX.size(); i++) {
282 if (cummulative) h+=histY[i];
312template <
class InputIt,
class OutputIt>
313inline void jkqtpstatKDE1DAutoranged(InputIt first, InputIt last, OutputIt KDEXOut, OutputIt KDEYOut,
double binWidth,
const std::function<
double(
double)>& kernel=std::function<
double(
double)>(&
jkqtpstatKernel1DGaussian),
double bandwidth=1.0,
bool cummulative=
false) {
314 double minV=0, maxV=0;
316 jkqtpstatMinMax<InputIt>(first, last, minV, maxV,
nullptr,
nullptr, &N);
318 std::vector<double> histX;
319 std::vector<double> histY;
321 const double binw=binWidth;
324 for (
double xi=minV; xi<=maxV; xi+=binw) {
328 if (histX.size()>0 && histX[histX.size()-1]<maxV) {
329 histX.push_back(maxV);
336 for (
size_t i=0; i<histX.size(); i++) {
338 if (cummulative) h+=histY[i];
367template <
class InputIt,
class BinsInputIt,
class OutputIt>
368inline void jkqtpstatKDE1D(InputIt first, InputIt last, BinsInputIt binsFirst, BinsInputIt binsLast, OutputIt KDEXOut, OutputIt KDEYOut,
const std::function<
double(
double)>& kernel=std::function<
double(
double)>(&
jkqtpstatKernel1DGaussian),
double bandwidth=1.0,
bool cummulative=
false) {
369 double minV=0, maxV=0;
371 jkqtpstatMinMax<InputIt>(first, last, minV, maxV,
nullptr,
nullptr, &N);
373 std::vector<double> histX;
374 std::vector<double> histY;
378 for (
auto it=binsFirst; it!=binsLast; ++it) {
381 std::sort(histX.begin(), histX.end());
384 for (
auto it=histX.begin(); it!=histX.end(); ++it) {
391 for (
size_t i=0; i<histX.size(); i++) {
393 if (cummulative) h+=histY[i];
423template <
class InputIt,
class OutputIt>
424inline void jkqtpstatKDE1D(InputIt first, InputIt last,
double binXLeft,
double binXDelta,
double binXRight, OutputIt KDEXOut, OutputIt KDEYOut,
const std::function<
double(
double)>& kernel=std::function<
double(
double)>(&
jkqtpstatKernel1DGaussian),
double bandwidth=1.0,
bool cummulative=
false) {
425 double minV=0, maxV=0;
427 jkqtpstatMinMax<InputIt>(first, last, minV, maxV,
nullptr,
nullptr, &N);
429 std::vector<double> histX;
430 std::vector<double> histY;
434 for (
double x=binXLeft; x<=binXRight; x+=binXDelta) {
442 for (
size_t i=0; i<histX.size(); i++) {
444 if (cummulative) h+=histY[i];
479template <
class InputItX,
class InputItY>
480inline double jkqtpstatEvaluateKernelSum2D(
double x,
double y, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY,
const std::function<
double(
double,
double)>& kernel,
double bandwidthX,
double bandwidthY) {
485 for (; (itX!=lastX)&&(itY!=lastY); ++itX, ++itY) {
489 const double vvx=(x-vx)/bandwidthX;
490 const double vvy=(y-vy)/bandwidthY;
491 res+=kernel(vvx,vvy);
495 if (cnt==0)
return 0.0;
496 return res/
static_cast<double>(cnt)/sqrt(bandwidthX*bandwidthY);
525template <
class InputItX,
class InputItY,
class FF>
526inline int jkqtpstatStoreKernelSum2D(FF&& fStore, InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY,
const std::function<
double(
double,
double)>& kernel,
double bandwidthX,
double bandwidthY) {
530 for (; (itX!=lastX)&&(itY!=lastY); ++itX, ++itY) {
534 fStore(vx,vy, kernel,bandwidthX, bandwidthY);
565template <
class InputItX,
class InputItY,
class OutputIt>
566inline void jkqtpstatKDE2D(InputItX firstX, InputItX lastX, InputItY firstY, InputItY lastY, OutputIt histogramImgOut,
double xmin,
double xmax,
double ymin,
double ymax,
size_t xbins,
size_t ybins,
const std::function<
double(
double,
double)>& kernel=std::function<
double(
double,
double)>(&
jkqtpstatKernel2DGaussian),
double bandwidthX=1.0,
double bandwidthY=1.0) {
568 const double binwx=fabs(xmax-xmin)/
static_cast<double>(xbins);
569 const double binwy=fabs(ymax-ymin)/
static_cast<double>(ybins);
572 auto itOut=histogramImgOut;
573 for (
size_t by=0; by<ybins; by++) {
574 for (
size_t bx=0; bx<xbins; bx++) {
581 const double cnt=
jkqtpstatStoreKernelSum2D([&](
double vx,
double vy,
const std::function<
double(
double,
double)>& kernel,
double bandwidthX,
double bandwidthY){
583 auto itOut=histogramImgOut;
586 for (
size_t by=0; by<ybins; by++) {
588 for (
size_t bx=0; bx<xbins; bx++) {
589 const double vvx=(x-vx)/bandwidthX;
590 const double vvy=(y-vy)/bandwidthY;
591 *itOut+=kernel(vvx,vvy)/sqrt(bandwidthX*bandwidthY);
598 }, firstX, lastX, firstY, lastY, kernel, bandwidthX, bandwidthY);
602 auto itOut=histogramImgOut;
603 for (
size_t by=0; by<ybins; by++) {
604 for (
size_t bx=0; bx<xbins; bx++) {
628template <
class InputIt>
632 return sigma/pow(
static_cast<double>(N), 1.0/(2.0+4.0));