【python因果推断库16】使用 PyMC 模型进行回归拐点设计

目录

例子 1 - 连续分段线性函数


import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import causalpy as cp
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
seed = 42
rng = np.random.default_rng(seed)

回归拐点设计应当通过分段连续函数来进行分析。具体来说:

  • 我们想要一个函数能够捕捉拐点左边和右边的数据。
  • 我们想要一个具有一个断点或拐点的分段函数。
  • 该函数应在拐点处连续。

一个这样的函数的例子是一个分段连续多项式:

\mu=\beta_0+\beta_1\cdot x+\beta_2\cdot x^2+\beta_3\cdot(x-k)\cdot t+\beta_4\cdot(x-k)^2\cdot t

其中:

- \beta 是未知参数,
- x 是连续变量,
- t 是一个指示变量,如果 x>=k 则取值为1,否则为 0
- k 是 x 在拐点处的值。

我们可以通过绘制一系列随机选择的 \beta 系数的函数图像来直观地展示这些函数的样子,但所有这些函数都具有相同的拐点在0.5。

def f(x, beta, kink):
    return (
        beta[0]
        + beta[1] * x
        + beta[2] * x**2
        + beta[3] * (x - kink) * (x >= kink)
        + beta[4] * (x - kink) ** 2 * (x >= kink)
    )


def gradient_change(beta, kink, epsilon=0.01):
    gradient_left = (f(kink, beta, kink) - f(kink - epsilon, beta, kink)) / epsilon
    gradient_right = (f(kink + epsilon, beta, kink) - f(kink, beta, kink)) / epsilon
    gradient_change = gradient_right - gradient_left
    return gradient_change


x = np.linspace(-1, 1, 1000)
kink = 0.5
cols = 5

fig, ax = plt.subplots(2, cols, sharex=True, sharey=True, figsize=(15, 5))

for col in range(cols):
    beta = rng.random(5)
    ax[0, col].plot(x, f(x, beta, kink), lw=3)
    ax[1, col].plot(x, np.gradient(f(x, beta, kink), x), lw=3)
    ax[0, col].set(title=f"Random  {col+1}")
    ax[1, col].set(xlabel="x")

ax[0, 0].set(ylabel="$y = f(x)$")
ax[1, 0].set(ylabel=r"$\frac{dy}{dx}$");

回归拐点分析的思想是拟合一个合适的函数到数据上,并估计在拐点处函数梯度是否有变化。

下面我们将生成一些数据集并演示如何进行回归拐点分析。我们将使用一个函数来生成具有所需特性的模拟数据集。

def generate_data(beta, kink, sigma=0.05, N=50):
    if beta is None:
        beta = rng.random(5)
    x = rng.uniform(-1, 1, N)
    y = f(x, beta, kink) + rng.normal(0, sigma, N)
    df = pd.DataFrame({"x": x, "y": y, "treated": x >= kink})
    return df

例子 1 - 连续分段线性函数

在这个例子中,我们将坚持使用一个简单的连续分段函数。

kink = 0.5
# linear function with gradient change of 2 at kink point
beta = [0, -1, 0, 2, 0]
sigma = 0.05
df = generate_data(beta, kink, sigma=sigma)

fig, ax = plt.subplots()
ax.scatter(df["x"], df["y"], alpha=0.5)
ax.axvline(kink, color="red", linestyle="--")
ax.set(title=f"Change in gradient at kink point: {gradient_change(beta, kink):.2f}");

我们可以使用常规的 cp.pymc_models.LinearRegression 模型,并通过巧妙地通过 formula 输入指定设计矩阵来强制执行连续分段特性。

在这个例子中,将公式设置为 "y ~ 1 + x + I((x-0.5)*treated)"(其中 0.5 是拐点)等同于以下模型:

\mu=\beta_0+\beta_1\cdot x+\beta_3\cdot(x-k)\cdot t 

拐点两侧的梯度变化是通过数值方法评估的。epsilon 参数决定了用于计算梯度变化的拐点两侧的距离。

result1 = cp.pymc_experiments.RegressionKink(
    df,
    formula=f"y ~ 1 + x + I((x-{kink})*treated)",
    model=cp.pymc_models.LinearRegression(sample_kwargs={"random_seed": seed}),
    kink_point=kink,
    epsilon=0.1,
)

fig, ax = result1.plot()

如果你想绘制推断出的梯度变化的后验分布,你可以按照以下方式进行。

ax = az.plot_posterior(result1.gradient_change, ref_val=gradient_change(beta, kink))
ax.set(title="Gradient change");

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/878832.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

Android中的冷启动,热启动和温启动

在App启动方式中分为三种:冷启动(cold start)、热启动(hot start)、温启动(warm start) 冷启动: 系统不存在App进程(App首次启动或者App被完全杀死)时启动A…

AcWing算法基础课-789数的范围-Java题解

大家好,我是何未来,本篇文章给大家讲解《AcWing算法基础课》789 题——数的范围。本文详细解析了一个基于二分查找的算法题,题目要求在有序数组中查找特定元素的首次和最后一次出现的位置。通过使用两个二分查找函数,程序能够高效…

Mysql InnoDB 存储引擎简介

InnoDB 存储引擎是 Mysql 的默认存储引擎,它是由 Innobase Oy 公司开发的 Mysql 为什么默认使用 InnoDB 存储引擎 InnoDB 是一款兼顾高可靠性和高性能的通用存储引擎 在 Mysql 5.5 版本之前,默认是使用 MyISAM 存储引擎,在 5.5 及其之后版…

车型展示+接驳体验!苏州金龙海格客车闪耀汉诺威商用车展

德国当地时间9月16日,IAA汉诺威商用车展媒体日活动在德国汉诺威展览中心开幕。该展会自1897年首次举办以来,已有超过一个世纪的历史,是全球历史最长、规模最大、最具影响力的专业商用车展之一,更是世界商用车行业技术创新和发展趋…

实战案例(5)防火墙通过跨三层MAC识别功能控制三层核心下面的终端

如果网关是在核心设备上面,还能用MAC地址进行控制吗? 办公区域的网段都在三层上面,防火墙还能基于MAC来控制吗? 采用正常配置模式的步骤与思路 (1)配置思路与上面一样 (2)与上面区…

STM32 如何生成随机数

目录 一、引言 二、STM32 随机数发生器概述 三、工作原理 1.噪声源 2.线性反馈移位寄存器(LFSR) 3.数据寄存器(RNG_DR) 4.监控和检测电路: 5.控制和状态寄存器 6.生成流程 四、使用方法 1.使能随机数发生器 …

C++笔记---二叉搜索树

1. 二叉搜索树的概念 二叉搜索树又称二叉排序树,它或者是一棵空树,或者是具有以下性质的二叉树: • 若它的左子树不为空,则左子树上所有结点的值都小于等于根结点的值。 • 若它的右子树不为空,则右子树上所有结点的值都大于等于…

数据结构—双向链表

结构 带头链表里的头结点&#xff0c;实际为“哨兵位”&#xff0c;哨兵位结点不存储任何有效元素&#xff0c;只是站在这里“放哨 的” 实现双向链表 List.h #pragma once#include<stdio.h> #include<stdlib.h> #include<assert.h> #include<stdbool…

“RISCV+AI”

概述 设计方案 主要有两种设计方案。 RISCV核ASIC RISCV核是标准的基于RISCV指令集的CPU设计&#xff0c;ASIC部分通常是基于RISCV自带的向量扩展指令集构建的向量处理器&#xff0c;或是自定义的矩阵计算单元。 根据CPUAI ASIC部件的接口可以分为紧耦合和松耦合的设计1。 …

基于python+django+vue的学生管理系统

作者&#xff1a;计算机学姐 开发技术&#xff1a;SpringBoot、SSM、Vue、MySQL、JSP、ElementUI、Python、小程序等&#xff0c;“文末源码”。 专栏推荐&#xff1a;前后端分离项目源码、SpringBoot项目源码、SSM项目源码 系统展示 【2025最新】基于协同过滤pythondjangovue…

【Python笔记】PyCharm大模型项目环境配置

一、PyCharm创建新项目 二、更新pip版本 ...>python.exe -m pip install --upgrade pip 三、生成所需requirements配置文件 ...>pip freeze > requirements.txt 四、安装所需组件requirements.txt ...>pip install -r requirements.txt

【Kubernetes】linux centos安装部署Kubernetes集群

【Kubernetes】centos安装Kubernetes集群 1、环境准备 系统centos7 配置yum源参考文章 Centos系统换yum源 yum -y update 步骤1-3是所有主机都要配置&#xff0c;主机名和hosts配置完后可以使用工具命令同步 1.1 主机 一主二从 主机名ipk8smaster192.168.59.148k8snode11…

Node.js 安装及项目实践

node.js安装 node安装&#xff0c;选择版本 一直next&#xff0c;或者自己修改路径&#xff0c;添加两个包 选择自己的安装的node的路径&#xff0c;cmd或者winr cmd 显示node与npm的版本号 node -vnpm -v可以跟着这个博客将node安装 2024最新版Node.js下载安装及环境配…

ZW3D二次开发_UI_非模板表单_设置表单显示位置

1.ZW3D弹出非模板表单时可以设置弹出位置&#xff08;居中、左下角、右上角等&#xff09; 2.假设已创建好非模板表单 3.在Form属性中添加form_pos属性 4.输入值 base,CTR,0.0 &#xff0c;如下图 也可以设置为其他值显示在不同的位置&#xff0c;如下 5.重新编译&#xff0c;…

新升级|优化航拍/倾斜模型好消息,支持处理多套贴图模型!

【天元轻量化软件】一直在不断地追求进步和完善&#xff0c;以满足更多用户的各种需求。 电脑登录天元官网免费体验&#xff1a;天元轻量化软件官网 本次我们对“智能PBR”功能进行了更新。更新后的“智能PBR”支持带多套贴图的模型进行使用。 本轮更新后&#xff0c;主要受益…

防火墙--NAT技术,基于源NAT,NAT服务器,双向NAT

文章目录 防火墙--NAT技术一、基于源NAT**方式**&#xff1a;NAT No-PATNAPT出接口地址方式Smart NAT三元组 NAT 二、基于服务器的NAT多出口场景下的NAT Server 三、双向NAT 防火墙–NAT技术 基于源NAT&#xff1a;用于将内部网络的私有IP地址转换为公共IP地址&#xff0c;以便…

51单片机应用开发---数码管的控制应用

实现目标 1、掌握数码管结构、驱动原理及应用&#xff1b; 2、掌握数码管段码表推导&#xff1b; 3、会编程让开发板8个数码管动态显示。 一、什么是数码管&#xff1f; 1.数码管定义 数码管&#xff0c;也称为LED数码管&#xff0c;基本单元是发光二极管(LED)。分为七段数…

【机器学习】--- 自监督学习

1. 引言 机器学习近年来的发展迅猛&#xff0c;许多领域都在不断产生新的突破。在监督学习和无监督学习之外&#xff0c;自监督学习&#xff08;Self-Supervised Learning, SSL&#xff09;作为一种新兴的学习范式&#xff0c;逐渐成为机器学习研究的热门话题之一。自监督学习…

某思CMS V10存在SQL注入漏洞

Fofa: product"魅思-视频管理系统" 框架:ThinkPHP 5,6 1 漏洞分析&复现 位于 /controller/Api.php 控制器中的getOrderStatus 方法POST传入&#xff0c;然后直接拼接了 orderSn 变量到 where 查询中&#xff0c;导致漏洞产生. /** * 查询订单支付状态 */ pub…

服务器——装新的CUDA版本的方法

服务器——装新的CUDA版本 一、进入 CUDA 版本列表二、根据自己服务器&#xff0c;选择对应的版本和配置三、使用管理员用户&#xff0c;运行下载和安装命令四、查看显卡驱动是否安装4.1 若安装了显卡驱动4.2 若显卡驱动没安装 参考文章 一、进入 CUDA 版本列表 CUDA Toolkit …