数据分析

一个关于风险评测计算的问题

前言

近几天，在做一个工作相关的内容，就是对企业投资者进行风险评测。我设计了若干个问题，并根据投资者的选项用于评分，但需要提前计算多少种选项的可能性，并计算分数的分布。并根据分数确定相应的分数区间对应的风险等级。

这听上去是个简单的问题，但实际上计算起来颇为复杂，想了不少办法才实现。这里也记录下来，并写成分享的文章可以让大家看一下，说不定可以借鉴。

本次计算主要是使用Python和R。

由于代码比较简单，就不用同步到Github，所有代码在文章里都能获取到。

问题

其实这次的问题很简单:

  1. 有6个风险评测的问题
  2. 每个问题有A、B、C三个选项
  3. 每个问题的分值比重不一样
  4. 每个选项的分数不一样
  5. 需要根据投资者的选择答案计算其最终得分
  6. 根据其得分给出风险等级

但我们要根据所有可能的选择结果确定分数范围，并根据分数范围确定其对应的风险等级。所以问题的计算逻辑就如下了：

  1. 获取所有的选择答案的组合结果
  2. 根据组合计算每一个结果的分数
  3. 根据所有分数确定分数区间和风险等级的对应关系
  4. 计算投资者的分数
  5. 确定投资者的风险等级

获取所有组合

我们假设有A、B、C三个选项，对应的分数分别为[25,50,75]:

[A, B, C] = [25, 50, 75]

总共有6道题，每道题有三个选项，所以答案的组合总共有$3^6=729$个，那么分数的可能性就有729个（实际上应该小于729，因为有不同答案但相同分数）。

我们不需要获取ABC的答案组合，只需要获得分数[25, 50, 75]的组合即可，所以我思考了不超过10s，就放弃了使用Excel，因为这么多年早忘了当初那丁点的VBA知识了，而这个要用Excel计算的话必须使用VBA，所以果断放弃。

那么再思考下要怎么计算，我们知道有6个答案，每个答案的分数是[25, 50, 75]中的一个，然后组成一个组合。所以其答案可能是:

[A, A, B, B, C, A] = [25, 25, 50, 50, 75, 25]

这样的组合总共有729个，所以要获取所有可能的结果，那么可以使用Python的itertools包中的product就是笛卡尔积来计算所有的可能组合：

import itertools

s = []
for i in itertools.product([25,50,75], repeat =6 ):
	s.extend([i])

s就是所有的答案分数的组合可能，我们查看下组合的数量：

print(len(s))
729
[Finished in 0.1s]

如何计算？

我们已经知道了，每个答案组合有对应的分数，以上面的组合为例：

[A, A, B, B, C, A] = [25, 25, 50, 50, 75, 25]

计算其分值为：250(随机的分数都能是这个值我服了)。

但我们不能这样简单的计算分值，因为每个问题的分值权重是不一样的，每个问题的重要程度是不一样的，我们预先定义每个问题权重为：

[1, 2, 3, 4, 5, 6] = [20%, 20%, 20%, 20%, 10%, 10%]

所以每个投资者的分数计算公式为：

$$Score = \sum_{i}^{n=6}score_i \times W_i$$

$score_{i}$即是答案的分值组合，$W_i$为每个答案分值的权重,重新计算的分数结果就是40分。

如果答案组合只有一两个，那么这样计算也就罢了，但我们有729个组合，不可能每个计算出来，至少我们的组合的获取就是要通过工具得出来。怎么计算所有729个组合分数呢，而且当选项扩展至4个甚至5个，问题扩展至数十个，那么答案组合就更多了，怎么办呢？

复习下矩阵计算

在6个问题3个答案的情况有729个组合，每个组合是有6个分数组成，每个分数有对应的权重，我们想到这样一个公式，$A_{m \times n}$ $\times$ $B_{n \times k}$ = $C_{m \times k}$:

$$
\begin{equation}
\begin{bmatrix}
1 & 2 \
3 & 4
\end{bmatrix}
\times
\begin{bmatrix}
1 \
2
\end{bmatrix}

\begin{bmatrix}
5 \
11
\end{bmatrix}
\end{equation}
$$

这是线性代数中所学习的矩阵乘法，那么我们就可以认为，729个组合每个组合为6个分数，所以可以看成一个$729 \times 6$的矩阵$A_{729 \times 6}$，那么权重则是一个$6 \times 1$的矩阵$B_{6 \times 1}$，所有两个矩阵相乘则生成一个矩阵$C_{729 \times 1}$。

那么计算的整个公式就是：

$$
\begin{equation}
\begin{bmatrix}
a_{11} & a_{12} & a_{13} & \cdots & a_{16} \
a_{21} & a_{22} & a_{23} & \cdots & a_{26} \
\vdots & \vdots & \vdots & \ddots & \vdots \
a_{i1} & a_{i2} & a_{i3} & \cdots & a_{ij}
\end{bmatrix}
\times
\begin{bmatrix}
b_{11} \
b_{21} \
\vdots \
b_{61}
\end{bmatrix}

\begin{bmatrix}
c_{11} \
c_{21} \
\vdots \
c_{i1}
\end{bmatrix}
\end{equation}
$$
特别的，$(i, j) = (729, 6)$。

Python计算逻辑和代码

我们考虑整个计算的顺序和逻辑如下：

  1. 获取所有分数组合
  2. 将分数组合生成为矩阵
  3. 构建权重的矩阵
  4. 分别两个矩阵相乘
  5. 得到所有分数
  6. 用R做分数的图形分布
  7. 根据分布确定不同风险等级的分位数

我们考虑每一个步骤的代码：

  1. 获取所有分数组合
import itertools

s = []
for i in itertools.product([25,50,75], repeat =6 ):
	s.extend([i])
  1. 将分数组合生成矩阵
import numpy as np

sm = np.mat(s)
  1. 构建权重的矩阵
w = np.mat([0.2,0.2,0.2,0.2,0.1,0.1])
  1. 矩阵相乘
score = np.matrix.tolist(sm*w.T)  % w.T指将矩阵转置为列矩阵

这段代码使用了np.matrix.tolist将计算后得到的矩阵直接转换为list列表，后面的第6步和第7步考虑在Python中调用R来实现。

Python调用R分析数据

数据的预处理

首先，需要使用python生成csv文件以供R使用，当然也可以在R中使用rpy2包调用R来分析。

  1. 使用python生成csv文件：
writer = csv.writer(file("score.csv",'wb'))
writer.writerow(['Value'])

for s in score:
	writer.writerow(s)
  1. 在R中导入数据：
rm(list = ls())
setwd("/Users/eggs/Library/Mobile Documents/com~apple~CloudDocs/文档/document")
score <- read.csv("score.csv",header = TRUE)

怎样分析数据？

其实分析的思路挺简单，总共有729个总体值，但实际上计算出来的Score值只有21个，所以我们要考虑每个Score值占总体组合的概率。一般我们考虑CDF（cumulative distribution function）就是累积分布函数：

$$
F(X) = F(x \leq X) =
\begin{cases}
\sum_{i=1}^{x \leq X}P_{x=i}， & \text{x是离散的} \[3ex]
\int_{-\infty}^{X} f(x),dx， & \text{x是连续的}
\end{cases}
$$

所以我们可以在R中画出Score的累积分布函数图。

画累积分布函数图

  1. 画出累积分布图
par(family='MicrosoftYaHei')
plot(ecdf(score$Value), do.points=FALSE, verticals=TRUE, main = 'CDF of score', xlab = 'score', ylab = 'CDF' )

画出来的图如下：
测试分数的累计分布图

  1. 用分位数工具取数

我们假设取[20%, 80%]分位数的值：

y <- quantile(score$Value,c(0.2,0.8))

 20%  80% 
42.5 57.5

那么，我们可以认为：

分数在[0,25)之间的为谨慎型，当然这一类别是理论值，不会通过问卷计算出来
分数在[25,42.5]之间的为保守型，占20%比例
分数在(42.5,57.5]之间的为谨慎型，占60%比例
分数在(57.5,75]之间的为积极型，占20%比例

结束语

在这里简单的做了一个关于风险评测的计算介绍，由于最近比较忙，所以内容写的也不怎么严谨。简单说就是思路清奇，排版华丽（你知道我说的是反话）。而且最近也没写什么，所以只是恰好碰到问题就随便写了篇，姑且看着吧。

python代码如下：

__author__ = "Eggs"
# coding: utf-8

import itertools
import numpy as np
import csv

s = []
for i in itertools.product([25,50,75], repeat = 6 ):
  s.extend([i])

sm = np.mat(s)
w = np.mat([0.2,0.2,0.2,0.2,0.1,0.1])

score = np.matrix.tolist(sm*w.T)

writer = csv.writer(file("score.csv",'wb'))
writer.writerow(['Value'])

for s in score:
  writer.writerow(s)

R的代码如下：

rm(list = ls())
setwd("/Users/eggs/Library/Mobile Documents/com~apple~CloudDocs/文档/document")
score <- read.csv("score.csv",header = TRUE)
a <- unique(score$Value)

par(family='MicrosoftYaHei')
plot(ecdf(score$Value), do.points=FALSE, verticals=TRUE, main = 'CDF of score', xlab = 'score', ylab = 'CDF' )
y <- quantile(score$Value,c(0.2,0.8))

length(score$Value)
length(a)
y

运行的结果如下：

> length(score$Value)
[1] 729
> length(a)
[1] 21
> y
 20%  80% 
42.5 57.5

一个测试：如果我们把题目答案扩展到4个，题目增加为10个，计算的结果如下：

> length(score$Value)
[1] 1048576
> length(a)
[1] 119
> y
 20%  80% 
43.5 56.5

可见总共的样本超过了100万，我们可以用图形观察如下：

CDF图形：
CDF

概率密度图形：
Denstiy

可见样本越多，图形越平滑，而且其分布越接近正态分布（极限中心定理）