求利用C++言进行的 lowess拟合 或者 loess拟合原代码

Posted 2023-04-13

前段时间用 R语言 对基因芯片数据进行Lowess标准化和loess标准化，但我提取数据都用的是C++做的，分析数据与提取数据用不同语言很麻烦，不知道那位大哥大姐，小弟小妹 有lowess或者loess拟合的C++原代码，若有能否给我参考一下下，在此谢过啦：）
loess拟合 又名 "局部加权回归法"

参考技术A 有 老 Fortran 语言源程序
进
http://www.stat.purdue.edu/~wsc/localfitsoft.html
取
LOWESS base software from netlib.
LOESS base software from netlib.

老 Fortran 翻译成 C 不难。 参考技术B ANN: LOWESS fitting (python wrapper)
Istvan Albert ialbert at mailblocks.com
Fri Sep 10 20:53:59 CEST 2004

Previous message: Solved: Dynamic Linking Problems with Intel Compiler
Next message: python to C ?
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

--------------------------------------------------------------------------------

Hello folks,

I have recently needed to use lowess fitting in a project. I ended
up using the C code from R by wrapping it as a python module.

If you need such things have a look here:

http://www.personal.psu.edu/staff/i/u/iua1/python/dist/

Windows binaries for python 2.3 are included in the library.

In the meantime I found out that biopython has a lowess call,
but that depends on Numeric. (It also does not seem to work
properly. As far as I could tell, it was not returning correct
values).

cheers,

Istvan.

-------------- README ------------------

Lowess fitting for python, version 0.5
======================================

ABOUT
-----

Lowess stands for "locally weighted regression". For more
information see:

http://www.itl.nist.gov/div898/handbook/pmd/section1/pmd144.htm

This module is a python wrapper around the lowess C code used
in the R project. It exposes a single function used as:

import lowess
ys = lowess.fit(x=X, y=Y, F=F, NSTEPS=NSTEPS, DELTA=DELTA)

that computes the smooth scatterplot of Y against X
using robust locally weighted regression. The function
returns the fitted values, computed at each of the
values of the horizontal axis in X. Parameter description:

X = abscissas of the points on the scatterplot;
VERY IMPORTANT: the values in X must be ordered
from smallest to largest!

Y = ordinates of the points on the scatterplot.

F = specifies the amount of smoothing; F is
the fraction of points used to compute each
fitted value; as F increases the smoothed values
become smoother; choosing F in the range .2 to
idea which value to use, try F = .5.

NSTEPS = the number of iterations in the robust
fit; if NSTEPS = 0, the nonrobust fit is
returned; setting NSTEPS equal to 2 should serve
most purposes.

DELTA = nonnegative parameter which may be used
to save computations. If let to the default value of
'None' DELTA will be esstmated as (max(X)-min(X))/50

For more information on the parameters see the
lowess_readme.txt file.

INSTALLATION
------------

On unix platforms SWIG is required. Type

python setup.py build

then

python setup.py install

On windows platforms copy the two files located in
the win32 directory to your python path, usually to the
c:\Python23\Lib\site-packages directory.

TESTING
-------

Run the runme.py file.

LICENSE
-------本回答被提问者采纳

最佳拟合线与R中的阈值

我正在对一些可以从RDD中受益的数据进行回归。因此，我想在x轴上的阈值0.5以上和以下显示最佳拟合/回归线。我正在努力做到这一点。我已经尝试过clip(x1,x2,y1,y2)命令，但它仍然在整个情节中绘制线条。我还尝试使用子集绘制回归线> / <0.5，这也给出了整个图的一条线。

用lowess线路可能会更好吗？这对我来说真的是一个未知的R领域，所以我真的不确定如何继续。

答案

没有示例数据集，很难说哪种方法最适合你，但你可以考虑geom_smooth中的ggplot。

library(ggplot2)

# simulated data
set.seed(123)
beta_low <- -3; beta_high <- 2
cut_off <- 0.5 
x = seq(0,1, by = 0.01)
y <- ifelse(x < cut_off, x*beta_low, x*beta_high) + rnorm(length(x),     
                                                    mean = 0, sd = 2)

# add a new variable, say "group", which indicates whether you're before the 
# cut-off or after
df <- data.frame(x, y, group = ifelse(x < cut_off, "before_cut", 
"after_cut"))

# geom_smooth in combination with the group argument in aes() will make sure 
# that lines are only shown for the specified regions >/< 0.5
ggplot(df, aes(x, y, group = group)) +
geom_point() + 
geom_smooth(method = "lm", fill = NA, fullrange = FALSE)

或者，base R解决方案：

part1 <- lm(y ~ x, subset=(x<=cut_off))
part2 <- lm(y ~ x, subset=(x>cut_off))
plot(x,y)
segments(min(x), part1$fitted.values[1],                             
         cut_off, rev(part1$fitted.values)[1])
segments(cut_off, part2$fitted.values[1],                             
         max(x), rev(part2$fitted.values)[1])